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Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

THE  ELECTRON  LIQUID  AT  ANY  DEGENERACY 
Brian  Gregory  Wilson 
May  1987 

Chairman:    Charles  F.  Hooper,  Jr. 
Major  Department:  Physics 

The  dielectric  formulation  of  the  many-body  problem  is  applied 
to  the  study  of  the  static  correlations  in  electron  liquids  at 
non-zero  temperatures.    The  intermediate  coupling  effects  arising  from 
the  exchange  and  Coulomb  correlations  are  taken  into  account  through  a 
local-field  correction  factor  obtained  from  quantal  extensions  of 
classical  fluid  integral  equations.    This  results  in  a  unified  theory 
which  can  describe  the  thermodynamic  and  dielectric  properties  of  the 
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electron  liquid  covering  the  limits  of  the  zero  temperature  jellium  to 
the  classical  one-component  plasma. 

Detailed  numerical  calculations  are  provided  for  the 
self-consistent  set  of  equations  describing  the  static  structure 
factor.    The  variation  of  interaction  energy  per  particle  versus 
temperature  at  constant  density  is  shown  to  exhibit  a  peaked  structure 
near  absolute  zero.    This  feature  is  not  observed  using  effective 
potential  approaches  in  classical  fluid  equations.    A  comparison  of 
zero  temperature  results  with  those  obtained  from  the  Greens  function 
Monte-Carlo  method  is  good;  however,  agreement  worsens  when  dynamical 
effects  of  the  local-field  correction  factor  are  included  using  the 
Mori  continued  fraction  method. 
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CHAPTER  I 
INTRODUCTION 


The  one  component  plasma  (OCP)  is  a  system  consisting  of  a  single 
species  of  charged  point  particles  together  with  a  uniform  oppositely 
charged  rigid  background  to  ensure  charge  neutrality.    It  was  first 
introduced  by  Sommerfeld  as  an  approximation  to  metals  when  band 
structure  effects  are  negligible.    By  coupling  the  standard  adiabatic 
(Born-Oppenheimer)  approximation  with  an  assumption  of  a  weak  electron- 
ion  pseudopotential  interaction,  it  can  be  shown^"^  that  the 
thermodynamics  of  a  simple  metal  is  the  sum  of  that  due  to  a  one 
component  electron  plasma  and  that  due  to  a  system  of  classical  ions 
interacting  via  a  state-dependent  pair  potential.    Viewed  as  a  model 
fluid,  the  OCP  exhibits  the  characteristics  which  distinguish  real 
coulomb  systems,  such  as  plasmas  and  electrolyte  solutions,  from 
neutral  fluids.    These  include  the  phenomena  of  screening  and  plasma 
oscillations,  which  arise  from  the  long-range  nature  of  the  coulomb 
interaction.    Thus  the  OCP  is  a  standard  approximation  for  both 
astrophysical  and  laboratory  fusion  plasmas,^  where  the  electrons  are 
treated  as  a  polarizable  background  for  the  gas  of  ions.  Further 
applications  lie  within  the  local  density  approximation  of  the  density 
functional  theory  of  electronic  structure,  where  the  excess  chemical 
potential  of  a  one-component  electron  plasma  plays  the  role  of  an 
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exchange-correlation  potential  in  the  one-electron  Schroedinger 
equation. 

It  is  thus  not  surprising  that  there  has  been  a  large  amount  of 
work  expended  in  obtaining  accurate  theoretical  descriptions  of  both 
the  zero  temperature  fermion  OCP  (degenerate  electron  liquid)  and  the 
high  temperature  classical  OCP  (electron  gas).^"^   Active  research  is 
ongoing  in  three  main  areas.    First,  in  the  course  of  study  of  linear 
dispersion  of  plasmons  in  the  simple  organic  polymers  polyacetyl ene 
and  polydiacetylene,  it  has  proved  important  to  contrast  the  properties 
of  the  exchange  hole,  and  its    screening,  in  quasi-one-dimensional 
solids  with  those  in  isotropic  three-dimensional  jellium  (classical 
OCP).^'^   Secondly,  one  has  been  able  to  determine  experimentally, 
using  inelastic  X-ray  scattering,  the  dynamical  structure  factor 
S(q,  w)  for  the  electrons  in  a  few  metals. ^^"^^    It  is  of  particular 
interest  that  the  experimental  S(q,  w)  show  a  tendency  of  two-peak 
structure  for  q  within  the  particle-hole  continuum,  a  feature  believed 
to  be  a  property  of  the  uniform  electron  liquid  not  presently  accounted 
for  accurately. ^^"^^    Lastly  renewed  attempts  are  being  made  to 
quantify  the  intermediate  degeneracy  regime  of  the  fermion  plasma. 

There  are,  in  fact,  many  physical  systems  whose  electrons  exist  in 
partially  degenerate  states  and  are  inadequately  described  by  either 
zero  temperature  or  classical  formalisms.    For  example,  the  finite 
temperature  thermodynamic  and  dielectric  properties  of  the  electron 
liquid  are  needed  for  a  proper  treatment  of  the  equation  of  state  of 
high  temperature  liquid  metals  forund  in  shock  wave  experiments^^ 
and  experiments  on  liquid  metals  expanded  to  near  their  critical 
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points,         of  high  density  inertially  confined  fusion  plasmas,  and 
of  the  finite  temperature  exchange-correlation  potentials  used  in 
density  functional  calculations  of  atomic  properties  at  high  pressure 
and  temperature.      Semidegenerate  electron  liquids  are  also  found 
in  semiconductors  where  the  low  densities  involved  lead  to 
correspondingly  low  Fermi  temperatures.    The  finite  temperature 
equation  of  state  of  the  electron  liquid  is  also  needed  for  a 
quantitive  explanation  of  the  miscibility  gap  in  solutions  of  alkali 
metals  in  their  alkal i-hal ide  melts. 

The  problem  of  the  OCP  at  intermediate  degeneracy  has  only 
recently  received  attention  because  of  complications  not  present  in 
the  extremes  of  a  classical  or  a  ground  state  quantum  description. 
For  example  at  zero  temperature  diagrammatic,  self-consistent 
distribution  function  methods^^"^^  and  pseudo  integral  equation 
methods, ^^"^'^  together  with  recent  Quantum  Monte  Carlo^^"^° 
calculations,  have  combined  to  present  an  accurate  description  of  the 
paramagnetic  electron  1  iquid.'*^"'*^    But  at  finite  temperature  we  no 
longer  deal  with  a  unique  ground  state  (invalidating  the  Monte  Carlo 
approach)  nor  may  we  utilize  the  ground  state  energy  variational 
principle  (invalidating  Fermi  Hypernetted  Chain  approaches).* 
Diagrammatic  approaches  may  be  generalized  to  finite  temperatures^^ 


Fermi  Hypernetted  Chain  calculations  of  quantum  fluids^  are 
actually  variational  calculations  of  the  ground  state  using  an 
adjustable  Slater-Jastrow  type  wavefunction  and  should  not  be 
confused  with  the  classical  integral  equation  method  per  se. Rather 
the  name  stems  from  a  diagramatic  expansion  of  the  ground  state 
distribution  function  which  is  regrouped  into  a  structure  of  the 
same  form  as  the  classical  HNC  equation. 
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but  are  perturbation  expansions  reliable  only  in  the  regime  of  weak 
44 

coupling. 

On  the  other  hand  the  classical  limit  may  be  handled  by  a  variety 
of  techniques,  including  the  Monte  Carlo  method/^"^°  the  molecular 
dynamics  method, and  various  approximate  methods  involving 
integral  equations,  most  notable  amongst  these  the  Hypernetted-Chain 
equation^^"^^  and  modifications  thereof .^^"^^   Here  quantal 
effects  introduce  two  complications.    First  a  quantal  calculation  of 
the  pair  distribution  function  is  hampered  fundamentally  by  the  fact 
that  the  canonical  partition  function  can  no  longer  be  factored  into  a 
solely  configuration-integral  part  and  a  part  involving  solely 
momentum  coordinates.    Secondly  the  thermodynamic  description  will 
also  depend  on  the  momentum  distribution  (or  off  diagonal  single  body 
density  matrix)^^"^^  and  this  information  is  not  trivially  contained 
in  the  pair  distribution  function. 

Various  attempts  have  been  made  to  overcome  these  complications 

within  the  framework  of  the  successful  methods  available  for  classical 

fluids.    Wigner-Kirkwood  expansions^^"^^  may  be  used  for  those  cases 

where  quantum  corrections  to  the  classical  partition  function  are 

perturbative,  but  their  convergence  as  a  series  in  temperature  or 

Plancks  constant  is  slow  (if  existent)  and  is  restricted  to  short 

ranged  potentials;^^  it  is,  in  fact,  an  ill  defined  expansion  for 
68 

the  OCP.       Published  work  often  centers  around  diffraction  effects 
and  assumes  Boltzmann  statistics. 

Previous  calculations  of  the  OCP  at  intermediate  degeneracy  and 
coupling  have  been  of  several  types.    The  most  obvious  course  has  been 
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to  approximate  the  true  quantum  mechanical  Slater  sum  which  appears  in 

the  canonical  partition  function  with  the  corresponding  classical 

expression  involving  a  temperature  and  density  dependent  effective 

potential.    Pokrant^^  and  others^^"^^  have  developed  a  finite 

temperature  variational  principle  to  approximate  the  Slater  sum  as  a 

product  of  an  ideal  fermion  Slater  sum  and  that  of  pairwise  additive 

potentials.    The  energy  is  calculated  by  approximating  three-body 

distribution  functions  in  terms  of  pair  distribution  functions,  which 

are  obtained  from  the  HNC  equation  using  the  effective  potentials. 

Zero  temperature  correlation  energies  from  this  method  are,  however, 

in  substantial  disagreement  with  the  parameterized  fittings  of  the 

quantum  Monte-Carlo  results. '^^"^^ 

Several  other  groups  have  studied  the  thermodynamic^^"^^*^^"^^ 

and  dielectric^^'''^  properties  of  the  OCP  by  using  finite 

temperature  perturbation  theory  within  the  random  phase 

approximation.    In  particular  Perot  and  Dharma-Wardana^^  and  Kanhere 
22 

et  al.     have  presented  closed  form  parameterizations  of  the 

exchange-correlation  energy  and  chemical  potential  within  the  random 

phase  approximation  (RPA)  for  the  paramagnetic  and  spin  polarized 

cases.    However  the  RPA  is  a  weak  coupling  theory  that  results  in 

considerable  error  at  metallic  densities. 

20 

Dandrea  et  al .     attempt  to  overcome  the  limitations  of  the  RPA 
by  including  a  static  local  field  correction  factor  (LFCF)  within  the 
framework  of  the  dielectric  formulation.    They  use  an  approximate  form 
for  the  LFCF  which  has  two  adjustable  parameters.    The  first  is  fixed 
by  assuming  a  value  for  the  pair  distribution  function  at  the  radial 
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origin.    The  second  is  obtained  from  an  appropriate  interpolation 
between  known  Monte-Carlo  results  at  zero  temperature^^ ""^^  and  in  the 
classical  regime^'^"^^  for  the  bulk  modulus.    Tanaka  et  al.'^^ 
approximate  the  LFCF  at  all  temperatures  by  the  static  zero 
temperature  LFCF  of  Singwi-Tosi-Land-Sjolander  (STLS)  self-consistency 
scheme. ^^"^^ 

This  thesis  concerns  itself  with  a  unified  theory  capable  of 
describing  the  dynamic  density  response  function  of  the  OCP  throughout 
a  broad  range  of  densities  and  temperatures:    from  a  fully  degenerate 
quantum  plasma  at  absolute  zero  temperature  to  a  nondegenerate 
classical  gas.    It  differs  from  the  approach  of  Dandrea  and  Ashcroft 
or  Tanaka  et  al .  in  two  important  respects.    First,  that  the  theory  be 
particularly  applicable  to  electron  liquids  at  metallic  densities, 
which  are  not  weakly  coupled,  perturbative  local  field  corrections  to 
the  random  phase  approximation  are  replaced  altogether  by  a  self- 
consistent  field  determined  by  a  quantal  extension  of  classical 
integral  equations.    Formally  we  can  define  a  LFCF  that,  unlike 
Dandrea  et  al . ,  requires  no  known  results,  and  unlike  Tanaka  et  al., 
is  temperature  dependent;  indeed  it  may  be  viewed  as  a  generalization 
of  the  STLS  LFCF  of  zero  temperature  and  reduces  to  the  proper 
classical  result.    Secondly,  dynamical  dependencies  of  the  LFCF  can  be 
rigorously  included.    This  has  been  shown  to  be  a  mechanism  for 
modeling  the  double  peaked  shape  of  the  dynamic  structure  factor 
observed  in  metals. ^^"^^    It  is  also  an  essential  theoretical  feature- 
static  LFCFs  cannot  simultaneously  satisfy  the  compressibility  sum  rule 
and  the  third  frequency  moment  sum  rule  for  the  density  response. 


CHAPTER  II 

BACKGROUND  THEORY 

This  thesis  will  concern  itself  solely  with  the  OCP  as  a  Fermi 

liquid.    The  most  basic  characteristic  of  a  liquid  is  that  it 

possesses  short  range  order  as  opposed  to  the  long  range  periodicity 

of  a  crystalline  solid.    Since  the  structure  of  a  crystal  is 

determined  experimentally  by  observing  the  Bragg  reflection  of  X-rays, 

it  is  natural  to  seek  a  quantitative  description  of  the  liquid 

structure  via  the  intensity  of  X-ray,  thermal  neutron,  or  light 
3  80—81 

scattering.  '  Of  particular  interest  then  are  the  fluctuations 

of  space  and  time  dependent  densities  which  describe,  at  long 
wavelength,  the  cooperative  motion  of  many  particle  systems.    For  this 
purpose  define  the  particle  density  correlation  function  as 

S^pCr  -  ?'.  t)  =  <6n(r,  t)  5n(r.  o)>^^ 

where  6n(r,t)  is  the  excess  particle  density  operator 

6n(r,  t)  =  n(r,  t)  -  <6n(r,  t)> 

eq 

and  the  equilibrium  average  is  indicated  by  the  brackets  <.  > 

eq- 

Specifically,  even  though  <n(r,t)>gq  =  0,  there  will  be  spontaneous, 
usually  small,  fluctuations  on  a  local  scale.    The  density  correlation 
function  is  a  measure  of  these  fluctuations. 
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In  this  chapter  we  will  begin  our  investigation  into  the 
structure  of  quantum  fluids  by  reviewing  the  general  properties  of  the 
time  correlation  functions  related  to  the  particle  density  correlation 
function.    These  include  1)  various  symmetries  (for  example  time 
reversal),  2)  positivity  properties  which  are  related  to  the  dynamical 
stability  of  the  system,  3)  fluctuation  dissipation  theorems  which 
connect  spontaneous  fluctuations  and  energy  dissipation  in  a  thermally 
equilibriated  system,  4)  Kramers-Kronig  relations  which  express 
causality,  like  the  well  known  one  between  the  index  of  refraction  and 
the  absorption  coefficient,  and  5)  sum  rules  which  provide 
restrictions  that  any  approximate  microscopic  theory  must  fulfill. 

The  System 

Because  we  will  be  dealing  with  a  many  particle  system  from  a 
quantum  statistical  point  of  view,  we  will  employ  the  second 
quantization  representation  of  quantum  mechanics  even  though  we  will 
mostly  be  working  in  a  canonical  ensemble  where  the  number  of 
particles  is  fixed. 

For  our  equilibrium  "system"  we  shall  define  in  the  Schroedinger 
(time  dependent  state)  representation  the  unperturbed  Hamiltonian 

H  =  T  +  V 

in  terms  of  field  creation  and  annihilation  operators. The 
kinetic  energy  operator  is 
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T  =  J  dl  ^""(I)  {- 


}  ^(I) 


while  the  interaction  energy  operator  is  given  by 


V  =  ^  J  ;  dl  dll  f^il)  <t*(U)  { 


}  *(II)  f(I) 


Roman  numerals  were  used  as  shorthand  notation  for  labels;  that 
is  "I"  stands  for  the  position  coordinate  as  well  as  the  discreet  spin 
label  of  particle  1 . 

If  there  is  an  external  (possibly  time  dependent)  field  present, 
there  is  in  addition  a  perturbation  to  the  hamiltonian  of  the  form 

6H(t)  =  -  ;  drn(r)  6V(r,  t) 
where  we  introduced  the  electron  density  operator 

n(r)  =  I 


In  the  external  field  the  hamiltonian  is  explicitly  time 
dependent  and  given  by 

H(t)  =  H  +  6H(t) 

in  the  Schroedinger  representation  where  the  operator  n(r)  is  time 
independent.    The  time  dependance  is  carried  by  the  density  matrix,  or 
ensemble  operator  p(t)  which  describes  the  state  of  the  system  such 
that  the  average  of  n(r)  is 

<n(r,  t)>  =  tr  p(t)  n(f)  tr  p(t)  =  1 


spins 


Basics  of  Linear  Response 
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What  follows  is  entirely  parallel  to  traditional  derivations  of  time 
dependent  perturbation  theory  in  elementary  quantum  mechanics  texts. 
We  have  to  solve  the  time  evolution  equation  for  the  density  matrix 

ifi  ^  p(t)  =  [H(t),  p(t)]  =  [H,  p(t)]  +  [6H(t),  p(t)] 

(note  the  commutator  is  reversed  from  the  Heisenberg  equation  of 
motion  for  operators)  subject  to  the  initial  condition 

P(t  =  — )  =  [H.  p^q]  =  0 

The  initial  condition  expresses  the  fact  that  the  system  is  stationary 

before  the  external  field  is  turned  on— we  also  require  that  the 

external  field  decay  sufficiently  rapidly  as  time  reaches  infinity. 

For  manipulations  it  does  not  matter  what  p     is  but  since 

eq 

the  system  starts  from  thermal  equilibrium  we  take 
Peq  =  P     ^^'^  P 

If  we  assume  a  time  dependance  to  the  density  matrix  of  the  form 
p(t)  =  pgq  +  6p(t) 

then  the  solution  to  the  time  evolution  equation  to  first  order  in 
6p  and  linear  in  the  external  field  is 

5p(t)  =  ^  /    dx  e-^"^t  -  -^^/^  [6H(x).  p    ]  e-^^'H^t  - 
_oo  eq 

(This  is  easily  verified  by  differentiation  remembering  that  time 
appears  both  in  the  limits  and  in  the  integrand.)    From  this  equation 
the  induced  change  in  the  average  density 
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6<n(rt)>  =  tr  p(t)  n(r)  -  tr  p  n(r) 
is  given  by 

5<n(r,  t)>  =  ;     df  J  dr'  <1  [n(r  t),  n(r'  t')]>     6V(r'  f ) 

_00  °Q 

t 

=  zi  ;     dt'  ;  df'  x"(ft;  f'  t')  6V(r'  t')  (2. A) 

where  [A,  B]  =  AB  -  BA  is  the  commutator  bracket.  In  the  above 
equation  are  the  Heisenberg  operator  for  the  unperturbed  system 

n(?'  t)  =  e^"*/^  nCr)  e"^"*/^ 

and  we  have  introduced  the  dynamic  densitv  response  function  (please 
note  the  double  prime  superscript)  defined  by 

x"(ft;  r'V)  [n(ft),  n(f't')]>gq 

This  equation  is  the  fundamental  result  of  linear  response 
theory.    It  shows  that  the  density  response  to  a  small  external 
potential  is  the  averaged  commutator  rather  than  the  correlation 
function  S(r.t)  as  one  might  expect.    However  we  shall  see  that  the 
two  are  intimately  related. 

The  Dvnamic  Densitv  Response  Fiinrtinn 


In  equilibrium  the  system  is  space  translational ly  invariant.  We 
adopt  the  convention  that  forward  spatial  Fourier  transforms  have  a 
minus  one  signature  and  write 
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x"(k)  =  J  d(?  -?)  e"^"^*^"^  "'''^  x"(r  - 

"  V  <  k  f^^t)'  n!(t')]> 
k  k 

Here  we  have  used  the  definition  of  the  spatial  transform  of  the  density 
operator 

n^(t)  =  ;  dr  e"^"^*^  n(?t) 
k 

which  can  be  expressed  in  terms  of  plane  wave  creation  and  annihilation 
operators  as 

njo)  =  i;  aj  ^a^  ^ 

from  which  it  follows  that 

n  ^  =  n^ 
-K  K 

Time  translation  invariance  of  the  dynamic  density  response 
function  is  easily  established  from  the  definition  of  a  Heisenberg 
operator  and  the  cyclic  nature  of  the  trace.    We  adopt  the  notation 
that  forward  temporal  Fourier  transforms  have  a  plus  one  signature: 

00 

x"(Kto)  =  ;   d(t  -  t')  e^"^*  -  ^'^  x"(K.  t  -  t') 

00 

It  is  easy  to  show  that  this  space  and  time  Fourier  transform  is 
real  (since  the  function  is  a  commutator  of  hermitian  operators),  an 
odd  function  of  frequency  (since  the  equilibrium  state  is  invariant 
under  time  reversal  and  parity),  and  depends  only  on  the  magnitude  of 
k  (spatially  isotropic).    It  can  also  be  shown  that 

«x"(Kco)  >  0 
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in  a  stable  system.      This  is  just  a  statement  that  a  dissipative 
many  body  system  takes  more  energy  out  of  the  external  field  than  it 
gives  back. 

In  addition  to  the  density-response  function  it  is  convenient  to 
introduce  certain  interrelated  functions.    The  complex  response 
function  (please  note  the  tilda)  is  defined  in  transform  space  as 

^(Kz)  =  f   ^  ^^^^ 
^  ir    CO  -  z 

-.00 

This  is  an  analytic  function  of  the  complex  variable  z  as  long  as 
imaginary  z  does  not  equal  zero  (on  the  real  axis  it  has  a  branch 
cut).    For  z  above  the  real  axis  it  reduces  to  the  temporal  Laplace 
transform  of  the  density  response: 

00 

x(Kz)  =  zi  ;   dt  e^^*  x"(Kt) 

0 

while  below  the  real  axis 

0         .  . 

=  (-  zi)  J     dt  e^^^  x"(Kt) 


This  follows  from  the  identity 

M  «i"t 
dco  e  

Ziri  co-z~  0  t<o 


"   f^     ^lut  izt 

r     dco_  e   e         t  >  0 

J  _  7     =  «      4.  ^  «        Inn  z  >  0 


—00 


0  t    >    0  T 

,'-,4.  Inn  z  <  0 

e"^^^    t  <  0 

which  can  be  proven  using  the  Cauchy  integral  formula  in  the  complex 
omega  plane.    (For  t>0  we  can  close  the  contour  on  top-exp(iwt) 
being  bounded-and  the  contour  is  in  the  positive  sense.    For  t<0  we 
must  close  the  contour  on  the  bottom,  the  contour  now  being  in  the 
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negative  sense.    We  then  just  sum  the  residues,  the  pole  being 
included  depending  on  the  sign  of  imaginary  z.) 

The  physical  response  may  be  defined  in  transform  space  as  the 
limit  of  the  complex  response  as  we  approach  the  real  frequency  axis 
from  above 


00 


x(Kco)  =  lim  x(K.  2  =  CO  +  ie)      =      limf     —    /"i'^"'^  , 
e-o  e-o    -co    ^   "  + 

Using  the  identity  (P  denoting  Cauchys  principal  value) 
=  P  1  ±  iir  5(x) 


the  physical  response  can  be  expressed  in  terms  of  its  real  and 
imaginary  parts  as 

x(ku)  =  x'(Kco)  +  ix"(Ku) 

where  the  imaginary  part  is  just  the  previously  defined  density 
response  and  the  real  part  (note  the  single  prime)  is 


x'(K(o)  =  PJ  dcoLx;^ 


_^    ir    u    -  0) 


an  even  function  of  frequency,  is  a  manifestation  of  the  Kramers- 
Kronig  relations. 

The  utility  of  the  physical  response  function  is  that  it  links 
the  induced  change  in  the  time-space  Fourier  component  of  the  density 
with  the  time-space  Fourier  component  of  an  arbitrary  external 
potential 

5<n(Kw)>  =  x(Ka))  6(Kco)  # 


-15- 


One  may  also  note  that  this  same  relation  may  be  easily  obtained  via 
eq.  2. A  using  the  definition  of  the  physical  response  in  real 
(untransformed)  space,  namely 

x(rt)  =  zi  e(t)  x"(r,t) 

The  response  to  a  static  external  potential  may  be  derived  from  linear 
response  theory  by  assuming  that  the  potential  is  adiabatical ly 
switched  on  starting  from  the  remote  past 

6V(rt)  =       «^<'^>  I  I  ° 

0  t  >  0 

At  t  =  0  the  external  potential  is  at  full  strength  and  we  find 

00 

6<n(ft  =  o)>  =  zi  J   dx  ;  dr'  x"(r  -  r',  x)  e""  6n(r') 

0 

or 

5<n(K,  t  =  o)>  =  x(K)  6V(K) 

where  the  static  response  is  given  by  zero  frequency  limit  of  the 
physical  response 

X(K)  =  x(K,  CO  =  0)  =  lim  x(K,  z  =  ie) 

e-O 

A  note  on  this  last  equality.    Formally  we  have 

lim  x(K.  z  =  ie)  =  lim  f   ^  =  (P  f   ^co  iHlMi^  ^  ^  ^  ^ 

e^O  e"*o  -00 

but  because  x"  is  a  real  odd  function  of  frequency  the  imaginary 
term  vanishes  and  we  can  drop  the  "principal  part"  condition. 
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Relation  to  Density  Correlation  Functions 


At  this  point  we  should  explicitly  confirm  the  statement  that  the 

density  response  is  intimately  related  to  the  density  correlation 

function,  the  latter  being  the  experimentally  pertinent  quantity.  Now 

because  the  Fourier  transforms  of  correlation  spectra  are  related  by 
83 

the  equation 

r  <AB(t)>_  e"^"*  dt  =  e"^^  f   <B(t)  A>     e"^"*  dt 

-00  _oo  °" 

it  follows  directly  that 

J     <|:^  [A.  B(t)]>     e-^"*  dt  =  j:^  (1  -  e"!^^)  J     <AB(t)>     e"^"^  dt 

and  so  in  particular  (using  the  fact  that  x"(K(o)  is  odd  in  frequency) 
x"(Kco)  =  ^  (1  -  e""^^)  SCKco) 

where  S(k,w)  is  the  time-space  Fourier  transform  of  the  time  dependent 
density  correlation  function.    We  note  that  S(r,t)  is  Fourier 
transformable  because  when  r  and/or  t  are  very  large,  then  n(r,t)  and 
n(0,0)  are  statistically  independent 


<n(r,  t)  n(Q,  0)>  „  =  <n(r,t)>  „  <n(0.  0)> 
*-M  "q  eq 

and  so  S(r,t)  as  defined  as  a  correlation  of  excess  density  vanishes. 

The  above  relation  between  the  density  response  and  density 
correlation  functions  is  an  explicit  example  of  what  is  called  the 
fluctuation  dissipation  theorem^^  which  relates  two  physically 
distinct  quantities  of  fundamental  experimental  significance:  1) 
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Spontaneous  fluctuations,  which  arise  even  in  the  absence  of  external 
forces  from  the  internal  motion  of  the  constituent  particles. 
Described  by  S(k,w)  these  fluctuations  give  rise  to  scattering  of 
neutrons  or  light.    2)  Dissipative  behavior— all  or  part  of  the  work 
done  by  external  stirring  forces—is  irreversibly  disseminated  into 
the  infinitely  many  degrees  of  freedom  of  thermal  systems. 

The  fluctuation-dissipation  theorem  shows  that  S(k,w)  is  not 
quite  symmetric  in  w,  being  a  little  stronger  at  positive  frequencies 
than  at  negative.    Indeed,  since  tox"(Ku)  is  always  even  in  w,  it 
i  s  general ly  true  that 

S(K,  -(o)  =  e"''^'^  S(K,  co) 

This  result  makes  sense  from  a  neutron  scattering  point  of  view. 
Positive  frequency  means  the  neutron  has  lost  energy  to  the  system  (by 
creating  an  excitation  of  energy  hw)  while  negative  frequency 
describes  a  process  in  which  the  neutron  has  picked  up  energy  from  the 
system  (by  destroying  an  excitation).    Of  course,  to  destroy  an 
excitation  you  must  first  have  one,  and  their  relative  abundance  is 
given  by  exp[-bhw].    This  dissymmetry  of  the  scattering  intensity, 
proportional  to  S(k,w)  is  only  pronounced  at  low  temperatures,  kt<hw 
—  it  is  absent  classically. 

7Q 

Further  Relations 

Besides  its  relation  to  the  time  dependent  density  correlation 
function  the  density  response  satisfies  certain  sum  rules,  usually 
presented  in  the  form  of  moments  of  the  density  response: 
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00 

«(o^»  E  ;  X"(K{o)}  do) 

_00 


The  j  =  0  or  zeroth  moment  follows  from  analyticity  of  x"(k,w). 
f    ^  =  X(K.  0)  =  X(K) 

_00 

All  odd  order  moments  vanish  because  x"  is  odd  in  frequency.  The 
remaining  higher  order  moments  can  be  derived  from  the  expression 

(^)''  x"(rt;  ?'t')  =  <^  [(i  1^)"^  n(r.  t)  ,  n(r't')]> 
taken  at  equal  times  t  =  t';  this  means  that 

r  x"(Kco)  =  J  d(?  -  f')  e^^*^''  -  '''^ 

_00 

<(j;)^''^  [[...[n(rt),  H]  ....  H] ,  n(r't)]> 

The  right  hand  side  contains  a  sequence  of  equal  time  commutators 
which  can  in  principle,  and  sometimes  in  fact,  be  exactly  calculated. 
For  example  the  second  order  moment  yields  the  Thomas-Rei che-Kuhn  or 
F-sum  rule  equivalent 

00  2 

r     do)   „ii/L^  \     e  K 
J     ^   wx  (Kco)  =  — — 

-CD  ^  "> 

For  the  electron  gas  the  fourth  moment  has  also  been  calculated. ^^"^^ 

The  sum  rules  are  often  used  to  provide  coefficients  for  an 
expansion  of  x(k,z)  for  large  z  from  the  definition 

vri^y'*  _  r"       x"(K(o)  _  ,«(/» 

-•'^Trco-z       -^     2     *   4~  +  — ^ 

_oo  2  z 
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From  its  derivation,  which  expands  (1  -  w/z)~   =  1  +  (w/z)  +  (w/z) 
+  ...  it  is  clear  this  expansion  can  only  be  asymptotic,  i.e.  valid 
only  when  Izl  is  large  compared  to  all  frequencies  in  the  system, 
which  means  all  frequencies  for  which  x"(k.,w)  is  not  effectively 
zero. 

An  interesting  feature  of  the  sum  rules  is  their  existence. 

There  is  no  reason  why  the  thermodynamic  average  of  the  multiple 

commutator  should  exist  to  all  orders,  and  indeed  although  the  sixth 

90 

moment  has  been  calculated  in  the  classical  fluid  limit     it  may 

91 

diverge  for  the  degenerate  electron  gas  in  certain  regions.  As 

the  sixth  moment  determines  the  w~^  term  in  an  inverse  frequency 

expansion,  its  divergence  would  herald  the  presence  of  a  term  of  the 
4  -11/2 

form  k  w         in  S(k,w)  at  long  wavelength  for  the  high  frequency 
behavior  of  the  spectrum.    Great  caution  is  thus  indicated  in  the  use 
of  moment  sum  rules  beyond  the  first  three  mentioned. 

Distribution  Functions 

It  is  natural  to  describe  the  structure  and  thermodynamic 
properties  of  liquids  in  equilibrium  by  employing  static  distribution 
functions.    Most  treatises  on  fluids  introduce  such  quantities  first 
and  later  generalize  to  include  time  dependance  as  a  prelude  to 
exploring  non-equilibrium  statistical  mechanics.    Here  we  have  taken  a 
reverse  path;  we  have  already  related  time  dependent  density 
correlation  and  density  response  functions  and  it  is  our  intent  now  to 
relate  these  with  the  static  distribution  functions  common  to  liquid 
theory. 
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The  most  basic  static  distribution  function  is  the  single 
particle  distribution  function.    If  we  define  the  number  density 
operator  in  N  particle  configuration  space  as 

n 

n(r)  =    I    6(f  -  r.) 

i  =  l  ^ 

(the  r.  are  quantum  operators)  then  the  single  particle  distribution 
is  given  by 

p(r)  =  <n(f)>  =  N  ^  J"  W(rp       ..  r^^)  dr2  ..  df^^ 
where  the  configurational  part  of  the  partition  function  is 
0  =  J  W(r^  . .  rj^)  df  ^  . .  drj^ 

and  the  configuration  probability  WCr^  ...  r^^)  is  given  by  the 
diagonal  slater  sum 


W(r^,  ...  r^Ir^,  ...  r^^) 


where,  in  general,  the  off  diagonal  Slater  sum  is  defined  as 


N(r^  -  r^  I  r|  ...r')  =  N  !  X^'^  Z  ?t(r^  ..  f^)  e"    °P  <i.Cr\ 


In  the  above  equation  the  index  "i"  denotes  any  complete  set  of 
properly  symmetrized  N-particle  basis  states  and  we  have  introduced 
the  DeBroglie  thermal  wavelength 

T  m 
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In  the  classical  limit  the  probability  function  "W"  is  just  the 
Boltzmann  factor  of  the  potential  energy  of  the  configuration;  at  zero 
temperature  it  consists  solely  of  the  ground  state  wavefunction. 

The  information  contained  in  the  single  particle  density 
distribution  is  minimal,  for  in  the  case  of  a  translational ly 
invariant  Hamiltonian  it  is  a  constant  equal  to  the  average  number 
density  of  the  system.    However  we  introduce  it  here  because  relevant 
thermodynamic  properties  we  will  consider  require  the  information  of 
the  off-diagonal  single  particle  distribution  function^"^"^^ 

p(r  I  r')  =  N  1  J  W(r^  r^...  r,^  I  r  j  r^  . .  r,^)  df^  ...  dr^ 

The  physical  significance  of  the  off-diagonal  distribution  function  is 
that  its  spactial  Fourier  transform  provides  us  with  the  momentum 
distribution  of  the  interacting  system. 

The  pair  distribution  function  is  defined  by 

g(r  ?')  =  N(N  -  1)  1  J  N  d?.  ..  r.  =  <Z  Z  Sir  -  f.)  6(f'  -  r,)> 

and  for  a  translational ly  invariant  system  is  of  the  form  g(r-r'). 
The  pair  distribution  function  in  turn  is  related  to  the  static 
structure  factor  through  the  relation 

S(K)  =  1  +  p  J  e^^*''  (g(r)  -  1)  dr 

Within  a  factor  of  density  the  static  structure  factor  is  simply  the 
time  independent  density  correlation  function. 

S(r  -?')=!  <(n(r)  -  <n(r)»  (n(r')  -  <n(r')»> 
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Note  that  if  we  had  defined  the  static  structure  factor  in  terms  of 
the  density  operator  instead  of  the  excess  density  we  would  replace 
g(r)  -  1  by  g(r)  in  the  formula  relating  the  two.    This  amounts  to 
including  the  experimentally  unmeasurable  "forward  scattering." 

The  definition  of  the  static  structure  factor  is  immediately 
generalized  to  form  the  dynamic  structure  factor,  essentially  the 
previously  introduced  time  dependent  density  correlation  function 
(again  within  a  constant  factor  of  density)  by  replacing  the  quantum 
operators  r.  with  their  time  dependent  equivalents  in  the  Heisenberg 
representation.    Bearing  in  mind  that  r.(t)  and  rj(t)  are 
operators  which  in  general  do  not  commute  (except  at  t  =  0),  we  easily 
obtain  several  important  properties  of  the  dynamic  structure  factor 
for  translational ly  invariant  systems: 

1)  S(-r,  t)  =  S*(r,t)  where  S*  is  the  complex  conjugate  of  S. 

2)  In  the  classical  limit  the  imaginary  part  of  S(r,  t)  vanishes. 

3)  S(K(o)  E^f   dt  e-^"^  ;  dr  e"^''*^  S(rt) 

_00 

is  real  in  both  the  classical  and  quantal  cases. 

4)  The  elastic  sum  rule 

00 

S(K)  =  J      S(Ka))  dto 

—CO 

is  an  immediate  consequence  of  the  definitions. 

All  the  relationships  amongst  the  various  functions  discussed  are 
best  conveyed  diagrammatical ly  (See  Fig.  2.1).    Note  that  in  the  upper 
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left  corner  we  have  introduced  a  new  function  C(k,w)  called  the  Kubo 
or  Relaxation  function,  about  which  we  will  have  much  to  say  later. 
For  now  it  simply  expresses  the  difference  between  the  (dynamic) 
physical  response  and  the  static  response  function. 

The  point  of  our  review  of  the  response,  correlation,  and 
distribution  functions,  so  concisely  summarized  in  Fig.  2.1,  is  that 
if  we  wish  to  generalize  the  classical  many-body  methods  of 
calculating  g(r)  into  the  quantal  regime,  there  are  two  possible 
relations  open  to  us.    First  there  is  s''^(k)  as  the  classical  limit 
of  $^""(1^)  (obviously)  or.  as 


P    6n(r  I  U) 
6  -  pU(f') 


=  -  ^  X(K) 


where  the  Iq  denotes  evaluation  at  constant  density/zero  external 
potential,  and  the  F  denotes  the  Fourier  transform  with  respect  to  the 
spactial  variable  r-r' . 

Let  us  contrast  the  two  alternatives.    The  quantum  mechanical 
static  structure  factor  S^'"(k)  is  essentially  statistical  in  nature; 
it  is  strongly  dependent  on  the  dynamics  of  the  particles,  depending 
on  S(k,w)  ~  Imag  x(Kw)  over  the  full  range  of  frequency.    On  the 
other  hand  x(K)  depends  on  the  single  zero  frequency  limit  of 
x(Ku),  and  is  essentially  mechanical,  arising  from  the 
Schroedinger  equation  by  our  linear  response  derivation,  which  makes 
but  implicit  reference  to  statistical  questions.    This  is  clearly 
indicated  by  the  F-sum  rule-which  holds  in  any  stationary 
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xCkco) 


X(kco)  =  x(K)  +  itoPCCku)  yiKta)  E  -  -  -  Im  x(Kco) 


ir  (0 


lim  (0  -  0  SCkco)  =  ^  (1  +  cot  h  ^)  Y(kco) 


x(K)  =  ;     y(Kw)  d(o 


X(K)  S(K)  =  J     S(Kto)  dco 


Figure  2.1 


Relations  between  the  linear  density-density  response  function  and  the 
structure  factor. 
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ensemble— expressing  nothing  but  the  conservation  of  particles  and  the 
elementary  commutator  of  position  and  momentum. 

This  dissertation  will  concern  itself  with  the  integral  equation 
methods  of  solving  the  many  body  problem.    This  route  seems 
particularly  promising  for  treating  quantal  fluids  in  that  for  the 
classical  limit  there  exists  a  wide  variety  of  integral  equations  for 
g(r)  which  have  been  investigated  and  give  results  in  good  agreement 
With  experimental  results.    '     The  purpose  here  is  to  extend  this 
procedure  into  the  quantal  regime.    A  fundamental  ansatz  of  this 
treatise  is  that  one  would  expect  the  more  accurate  integral  equation 
approximation  to  correspond  with  the  less  sensitive  function.  The 
extension  of  classical  integral  equations  into  the  quantal  regime  via 
the  density  response  function  is  presented  in  Chapter  II. 

With  this  fundamental  ansatz,  however,  there  is  a  fundamental 
drawback:    From  x(K)  one  obtains  an  integral  equation  not  for 
g^"^(r)  (which  is  used  in  the  computation  of  thermodynamic 
quantities)  but  rather  n(r/v),  that  is  the  distribution  function  in 
the  presence  of  an  external  potential  which  has  the  form  of  the 
interparticle  interaction.    While  it  is  true  that  in  the  classical 
limit  n(r/v)  =  n^qir) ,  this  is  not  true  in  the  quantal  case,  since 
we  are  neglecting  exchange  and  recoil  effects.    The  external  potential 
can  be  thought  of  as  arising  from  a  distinguishable  particle  of 
infinite  mass  and  therefore  fixed  in  position.    An  extreme  example 
occurs  for  the  case  of  very  weak  interactions  between  fermions.  In 
such  a  case  n(r/v)/nQ  is  approximately  unity  everywhere,  whereas 
g(r)  has  approximately  the  ideal  fermion  value  of  1/2  at  the  origin. 
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Thus  our  task  is  twofold:    first  to  generate  an  accurate  quantal 
integral  equation  for  n(r/v)/nQ  and  second  to  generate  an 
approximate  g(r)  from  knowledge  of  n(r/v)/n^.    In  essence  our  second 
task  is  to  go  from  the  bottom  left  of  Fig.  2.1  to  the  bottom  right 
along  the  paths  drawn  in  this  figure.    We  cannot  go  along  the 
diagonal,  as  this  path  is  one   directional  (we  lose  information  when 
we  integrate  y(k,vi)  over  w  in  order  to  obtain  x(k)).    We  must 
first  traverse  over  the  upper  left  path  of  Fig.  2.1;  we  must  obtain  an 
accurate  approximation  to  the  Relaxation  function.    This  is  taken  up 
in  subsequent  chapters. 


CHAPTER  III 
QUANTAL  FUNCTIONAL  EXPANSIONS 


In  this  chapter  we  will  establish  integral  equations     for  the 
quantity  n(r\y)/n^  (the  density  distribution  in  the  presence  of  an 
external  potential  of  the  form  of  the  interparticle  potential)  through 
the  functional  expansion  method. ^'^"^^    Because  in  the  classical 
limit  these  equations  reduce  to  the  familiar  Hypernetted-Chai n  (HNC) 
or  Percus-Yevick  (PY)  equations  for  the  pair  distribution  function,  we 
will  review  the  classical  derivation  with  the  aim  of  establishing 
corresponding  quantal  generalizations. 


The  OOZ  Relation 

In  classical  fluids  the  Ornstein-Zernicke  (OZ)  relation  plays  an 
important  role  in  treating  integral  equations  for  g(r).  Thus  in  this 
section  we  derive  a  quantal  extension  of  the  OZ  relation  (QOZ). 

We  start  from  the  known  functional  dependance  of  the  classical 
canonical  partition  function  on  pairwise  additive  potentials.  This 
allows  us  to  write  the  functional  derivative  of  the  density 
distribution  as^^ 

=  n^  6(?  -  ?■)  +  n^^  ^g(,^  _  r'l)  -  1}  =  F"^  {n^  S(K)} 

0 


6n(r  I  U) 
6  -  3U(r') 
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The      denotes  evaluation  when  the  arbitrary  static  external 
potential  U(r)  is  turned  off,  that  is  at  n(r)  =  n    (n  ;  the 

0  0 

average  density).    The  (negative  power)  script  F  denotes  the  (inverse) 
Fourier  transform  on  the  spatial  variable  r  -  r',  and  S(k)  is  the 
previously  defined  static  structure  factor 

S(K)  =  1  +  n^  ;  dr  e"^'^*'^  {g(r)  -  1} 
Now  consider  the  functional 


S  -  BU(r') 
6  N(f") 


This  is  a  vacuous  quantity  unless  such  a  functional  of  density  can 
actually  be  constructed.^^    If  such  a  creature  does  not  exist,  we 
define  it  through  the  chain  rule  of  functional  calculus: 


6  n(r") 


Sn(r") 
6  -  3U(r') 


=  6(r  -  r') 


It  follows  that  the  inverse  functional  derivative  is 


6  -  BU(r') 
6  n(r") 


1 


K    ^n    S(kT>  =  r         -        -  Cdr'  -  ?"l) 

0  0 


where  we  have  split  up  the  expression  for  the  derivative  into  two 
pieces  (this  should  be  considered  an  arbitrary  process)  with  the  first 
piece  chosen  to  be  the  non-interacting  (g(r)  =  1)  result  and  the 
second  piece  (c(r)),  called  the  direct  correlation  function, 
representing  the  deficit  from  from  the  non-interacting  result.  The 
functional  chain  rule  as  expressed  in  terms  of  c(r)  takes  the  form 


-29- 


h(r)  =  c(r)  +      J  c(lr  -  r' I)  h(r')  dr' 

This  is  the  classical  Ornstein-Zernicke  (OZ)  relation,  where  h(r)  = 
g(r)  -  1  is  called  the  pair  correlation  function. 

Furthermore,  in  the  classical  limit  where  there  are  no  exchange 
or  recoil  effects  (the  momentum  integrations  factor  from  the 
configurational  piece  of  the  partition  function)  one  can  rigorously 
derive  the  bootstrap  relation^^ 

n(r  I  U  =  v)  =  n^  g(r) 

where  the  external  potential  U(r)  has  the  form  of  the  interparticle 

interaction  V(r).  This  bootstrap  relation  allows  us  to  reformulate 

the  OZ  relation  in  terms  of  the  inhomogeneous  density  distribution  as 
fol lows: 


6  n(r') 


n(r'  I  U  =  v)  -  n^   =  c(r)  (3.1) 


0 

A  quantal  extension  of  the  OZ  relation  is  straightforward;^^  we  will 
follow  many  of  the  above  steps. 

First,  from  linear  response  theory,  we  saw  that  in  the  quantal 


case 


Sn(r  I  U) 
6  -  3U(r) 


^  =  I^K^  {%  X(K)} 


where  x(k)  differs  by  constant  factors  from  the  static  response 
function  previously  defined  in  Chapter  II.    Defined  here  it  is 
conveniently  dimensionl ess ;  note  that  only  the  product  (5x(k)  remains 
finite  as  temperature  goes  to  zero. 
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Again  the  chain  rule  of  functional  differentiation  requires  that 


6  -  BU(r') 
6n(r"  I  U) 


=  ^K^  ^=k^K    ^J(K)>  -  -  ^"'>  (3.2) 


but  here  we  have  decomposed  the  inverse  derivative  into  two  new  (and 
as  yet  unspecified)  functions.    These  are  determined  as  follows. 

As  we  wish  c(r)  to  reduce  to  the  classical  direct  correlation 
function  in  the  proper  limit,  we  need  only  require  J(k)  =  1  in  the 
classical  limit.    We  furthermore  make  the  ansatz  that 


J  ^  6n(r') 


n(r'  I  U  =  v)  -  n^   =  c(r)  (3.3) 


also  holds  for  the  quantal  case. 

We  have  seen  that  this  equation  is  strictly  true  in  the  classical 
case;  by  combining  the  two  above  equations  we  constrain  the  function 
c(r),  and  so  simultaneously  J(k),  by  the  relation 

F,{n(r  I  U  .  V)  -  n^)  =  xili  .  ,  ^^  ^^ 

In  hindsight  we  now  see  the  motivation  of  requiring  eq.  3.3  be 
fulfilled:    as  J(k)  tends  to  unity  in  the  classical  limit  we  have 
merely  imposed  the  condition  that  we  recover  the  bootstrap  relation  in 
the  classical  limit. 

The  above  equation  also  determines  J(k)  (and  so  c(r))  in  the 
general  quantum  limit,  for  if  the  system  described  by  n(r/v)  and 
x(k)  were  that  of  non-interacting  fermions,  then  in  the  limit  that 
U(r)  vanishes  we  see  that  J(k)  must  equal         .    This  is  an 
appealing  result,  for  it  amounts  to  defining  the  quantal  direct 
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correlation  function  by  separating  out  the  non-interacting  result, 
just  as  we  did  in  the  classical  case.    Henceforth  we  will  implicitly 
use  this  result,  which  will  be  referred  to  as  the  quantal  bootstrap 
relation  (QBS). 

As  a  consequence  of  this  now  specified  decomposition  of  the 
inverse  functional  derivative  eq.  3.2,  the  functional  derivative  chain 
rule  expressed  in  terms  of  the  quantal  direct  correlation  function 
becomes 


^  -  1)  =  M{c(?)}  +  n   ;  dr'(M{c(lr  -  r'l)})  (^i^^^  _  i) 
0  " 

where  M{...}  is  an  operator  defined  by 
F|^{Mf(r)}  =  Xq(K)  F^{f(?)} 

This  is  the  quantal  generalization  of  the  OZ  relation  (denoted  QOZ  for 
short).    Note  that  this  equation  can  be  cast  into  the  classical  form 
by  defining  analogue  functions: 

gan(-)  =  nirj_vi  ^an^-^  ^  ^an^-^  _  ^ 

0 

c^"(r)  E  M{c(r)} 

The  classical  equivalents  of  these  functions  all  have  "physical" 
interpretations  that  are  only  loosely  identified  with  the  analogue. 
(For  example  g^"  should  not  be  confused  with  the  true  or  quantum 
mechanical  pair  distribution  function.)    To  distinguish  classical 
analogues  we  shall  henceforth  employ  the  superscript  "an." 


-32- 


Classical  Generating  Functionals 


The  heart  of  the  functional  expansion  method  of  deriving  integral 
equations  lies  in  the  ability  to  construct  a  functional  W[n]  of  the 
density  distribution  n(rlu)  in  a  nonuniform  system  which  is  reasonably 
insensitive  to  the  function  n(r|u).^^    (For  other  points  of  view, 
see  Bakshi^°°  and  ChoquardJ^^)    If  such  is  the  case  we  can 
approximate  the  functional  W[n]  by  making  a  functional  "Taylor  series" 
expansion,  for  example  about  the  average  density  n^,  and  truncate 
after  the  first  (linear)  term 


W[n]  =  WCn^]  +  ;  dr'( 


6w 


6n(r') 


)  (n(r'  I  U)  -  n^) 


The  OZ  relation  then  forms  a  closed  set  of  two  equations  in  two 
unknown  functions.    As  an  illustration  we  will  consider  a  set  of 
functionals  and  the  integrals  equations  they  generate  in  classical 
statistical  mechanics. 

First  consider  the  functional 

W[n]  =  In         '  U> 
e-PU(r) 

The  denominator  forms  what  is  called  the  reference  distribution.  In 
order  to  make  the  functional  variational ly  insensitive,  it  is 
constructed  to  approximate  the  numerator  as  close  as  possible.  Here 
we  have  taken  the  classical  non-interacting  Boltzmann  factor.  The 
first  functional  derivative  is  easily  found  to  be 
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6w     ^  6(r  -  r')  _  6  -  SU(r) 
6n(r')       n(r  I  U)       6  n(r') 

and  from  the  definition  of  the  direct  correlation  function: 


6  -  m?) 
6n(r') 


6(r  -  r') 
=   ^  -  c(r  -  r' ) 

0  ° 


we  find 

Sw 


6n(r') 


=  c(r  -  r' ) 

0 


The  functional  expansion  of  W[n]  as  the  potential  U(r)  is  changed 
from  the  intitial  value  0  to  its  final  value 


In     "^^  '        =  ;  d?'(n(?'  I  V)  -  nj  (- 


6w 


) 


along  with  the  bootstrap  relation  n(r/v)/nQ  =  g(r)  and  the  OZ  relation 
then  combines  to  give  the  HNC  equation 

g(r)  =  e-  e^^"^^ 
where  yd)  =  h(r)  -  c(r)  is  called  the  nodal  function. 
The  PY  equation 

g(?)  =  e-  PV^^^  [1  .  Y(-r)] 
follows  from  the  same  steps  by  starting  with  the  functional 

n^  e- 

Because  the  HNC  functional  is  the  logarithm  of  the  PY  functional, 
it  is  the  more  insensitive  functional  of  the  two  from  variations.  One 
might  be  tempted  to  conclude  that  it  would  therefore  produce  the 
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superior  approximate  integral  equation.    However  the  validity  of  the 

HNC  equation  (or  any  other  derived  by  a  functional  expansion  method) 

depends  on  whether  the  cumulative  effects  of  the  higher  order 

derivatives  in  the  functional  expansion  are  negligible.    In  actuality 

their  effect  is  to  push  the  form  of  the  equation  towards  that  of  the 

PY  (note  how  the  form  of  the  PY  equation  is  the  small  y(r)  expansion 

of  the  exponential  y(r)  term  of  the  HNC).    This  is  evidenced  by  the 

fact  that  for  purely  repulsive  potentials  the  PY  and  HNC  equations  of 

state  bracket  the  "exact"  Monte-Carlo  or  molecular  dynamics  computer 

1 02 

simulation  results.        What  we  do  know  from  diagrammatic 
52—58 

analysis         is  that  the  exact  closure  equation  to  the  OZ  relation 
for  g(r)  is  of  the  form 

g(r)  =  e"  '^^^''^     Y^r)  +  B(r) 

where  the  unknown  "bridge  function"  B(r)  is  the  sum  of  bridge  diagrams 
of  two  point  functions.    The  evaluation  of  the  bridge  diagrams  in 
terms  of  higher  functional  deri vati ves^°^*^^  or  other  attempts  to 
account  for  bridge  effects  will  be  discussed  in  Chapter  V.    Here  our 
intent  is  to  generate  zeroth  order  integral  equations  from 
variationally  insensitive  functional,  thereby  minimizing  the  effects 
of  the  higher  order  terms  that  we  neglect  in  the  Taylor  functional 
expansion. 

It  is  well  known  that  the  PY  equation  has  no  solution  for  a  one 
component  plasma  system. Numerically  this  arises  from  the  fact 
that  in  the  PY  equation  for  g(r)  the  long  range  tail  of  the  potential 
is  not  compensated  by  that  of  the  nodal  function  yd),  as  it  is  in 
the  HNC  equation. This  is  symptomatic  of  a  physically  meaningful 
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inadequacy,  namely;  in  real  systems  the  particles  do  not  feel  the 
"bare"  long  ranged  potential  but  a  collectively  screened  or  short 
ranged  potential . 

An  improvement  on  the  PY  generating  functional  would  be  to 
replace  the  bare  potential  appearing  in  the  Boltzmann  factor  of  the 
reference  distribution  by  an  appropriately  screened  one.    For  example 

W[n]  =     "^'^  ' 

-  l5U„(r) 

"o  ^ 

where  the  Hartree  potential  is  defined  as 

3U^(r)  =  3U(r)  +  n^  /  dr'  (3U(r  -  ?')  (g(r')  -  1) 

As  the  Boltzmann  factor  now  more  closely  describes  the  physical 
distribution  n(r/v)  this  functional  is  variational ly  less  sensitive 
than  the  usual  PY 

_6w_  ^     Sir  -  ?■)    _     n(r  I  U)     ^  "  ^^h^"^^ 
6n(?')  -  3U.(r)  -  3Um(?)  6n(r') 

Using  the  bootstrap  relation  in  the  definition  of  the  Hartree 
potential  we  see  that 

LlJMIl  ^  s  -  m?)  ^  ^  ^5  (S  -  m?  - 

Sn(r')  6n(r')  6n(r') 

(n(S)  -  n^)  -  ;  dS  fJU(?  -  5)  6(?'  -  s) 

which  when  evaluated  in  the  limit  v  =  0(n(r|v)  =  n^)  and  using  the 
definition  of  the  direct  correlation  function  (eq.  3.4)  yields 
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6  -  3U^(r) 


=  ^^'^n  "^'^  -  c(r  -  r')  -  3U(r  -  ?') 
0  ° 


6n(r') 
and  so 

=  c(r  -  ?')  +  3U(r  -  ?') 


6n(r')  Q 

The  linear  expansion  of  the  functional  W[n]  thus  generates  the  equation 
"^^  '         =  1  +  ;  d?'  (n(?'  I  U)  -  n^)  (c(?  -  ?')  +  mr  -  ?')) 

which  can  be  re-expressed  with  the  help  of  the  bootstrap  equation  and 
the  OZ  relation  as 

-  BV.(?) 

g(r)  =  e  {1  +  pV^C?)  +  y(?)  _  pvCr)} 

Note  that  in  this  equation  the  long  range  part  of  the  potential 
is  cancelled  by  the  long  range  part  of  the  nodal  function,  and  all 
other  functions  appearing  in  this  equation  are  screened.  This 
equation,  which  we  shall  term  the  PYH  equation  (for  Percus-Yevi ck- 
Hartree)  is  equally  applicable  to  long  and  short  range  potentials, 
unlike  the  PY  equation. 

If  we  try  to  improve  the  HNC  equation  by  introducing  the  screened 
Hartree  potential  into  its  generating  functional 

W[n]  =  In  '  U) 

-  3U^(r) 

"o  ' 
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we  find  that  its  first  variation  is  the  same  as  in  the  above  PYH 
derivation  and  (following  many  of  the  same  steps)  the  functional 
expansion  leads  to  the  equation 


namely  the  original  HNC  equation.  In  other  words  the  HNC  equation  has 
the  striking  property  of  being  invariant  to  (Hartree)  screening  of  the 
reference  distribution  in  its  generating  functional. 

From  a  related  point  of  view  it  has  also  been  shown^°^  that  the 
HNC  equation  is  the  limit  in  a  series  of  integral  equations  that  can 
be  systematically  generated  starting  from  the  PY  equation,  among 
others. 

Finally  we  should  note  that  the  HNC  equation  can  be  derived  from 
an  even  more  interesting  point  of  view.    We  ignore  higher  order 
functional  derivatives  and  assert  that  the  best  functional  for  our 
closure  purpose  are  insensitive  to  variations,  and  this  property  shows 
up  in  the  magnitude  of  its  derivative.    This  can  be  minimized  by  using 
an  even  more  realistic  screened  potential  for  the  reference 
distribution.    In  fact  it  is  easy  to  show  that  if  we  modify  the 
Hartree  potential  by  replacing  the  convoluted  bare  potential  by  the 
direct  correlation  function 

W^ir)  =  (3U(r)  -  /  dr'  c(r  -  ?')  (g(?')  _  i) 
then  the  derivative  of  both  the  HNC  and  PY  like  functional 
H  =  In  -^U)  ,,        n(?  I  IJ) 


g(r)  =  e 


-  BV^(r)    Y(r)  +  3V^(r)  -  3V(r) 


n 


-  3Ujr) 


n 


0 


e 


-  3U  (?) 
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actually  vanishes.    They  in  turn  both  yield  the  same  equation 

-  3V  ^r) 
g(r)  =  e  ^ 

which  one  will  now  notice,  via  the  OZ  relation,  is  simply  the  HNC 
equation. 

Quantal  Generating  Functionals 

Corresponding  quantal  equations  can  be  obtained  from  the  above 
classical  generating  functionals  if  we  replace  the  reference 
distributions— the  Boltzmann  factor  n^expf-bU}  describing 
classical  non-interacting  particles  in  the  presence  of  an  external 
potential  U(r)— by  its  quantal  extension  n*(rlu)— the  density 
distribution  of  non-interacting  fermions  in  an  external  potential. 
(Henceforth  non-interacting  fermion  system  quantities  will  be 
distinguished  by  a  star.)    That  is  we  replace  the  reference 
distribution  by  a  single  particle  Hamiltonian  approximation.    As  in 
the  classical  presentation  of  the  above  section  the  functionals  are 
then  linearly  expanded  and  evaluated  in  the  limit  where  the  external 
potential  is  set  equal  to  the  interparticle  interaction  potential. 

We  have  previously  derived  a  quantal  extension  of  the  pair 
correlation  function 

V"o  ^("^^>  -rk)-^ 
from  the  ansatz  that  the  exact  classical  relation 


J  (6  -  SU(r) 
6n(r' ) 


)  (n(r'  I  U  =  V)  -  n^)  dr'  =  c(r) 
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remains  valid  in  the  quantal  case,  alternatively  expressed  by  the 
ansatz  that  the  classical  bootstrap  relation  is  extended  as 

By  making  use  of  these  relations  and  defining  the  analogue 
potential 

g-  3U^"(?)  ^  n*(?  I  U) 
"o 

a  closed  set  of  equations  is  formed  with  the  QOZ  relation.  Expressed 
in  terms  of  previously  introduced  analogue  functions,  these  integral 
equations  bear  a  close  correspondence  with  their  classical 
counterparts,  and  indeed  reduce  to  those  counterparts  in  the 
appropriate  limit. 

As  a  first  example  we  consider  the  functional 

W[n]  =  In    "^"^  I 

n*(r  I  U) 

Its  first  variation  is 

_-5w_  _  n*(r  I  U)  r6(?  -  r')        n(r  I  U)     Sn*(r  I  U), 
n(r  I  U)    ^n*(r  I  U)  "  ^^.^^  ,  ^^^2     6n(r")  > 

This  can  be  simplified  in  the  following  manner: 


6n(r')  *  5  -  3U(?")       ^     6n(r')  °      ^'''n^  xCQ) 


) 


=  1  -  FJc^"} 


where  it  should  be  remembered  that  the  analogue  pair  correlation 
function  differs  from  the  quantal ly  extended  "physical"  pair 
correlation  through  the  operator 
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-1 


M{f(r)}  =  F-'{x*(0)  F{f(r)}} 


From  the  above  considerations  we  see  that 


6n*(r  I  U) 
6n(r') 


=  6(r  -  r')  -  n^c^^Cr  -  r') 


and  so 


6w 


5n(r') 


=  c'"(?  -  r') 


By  Taylor  expanding  the  functional  W[n]  about  U(r)  =  0  to  v(r),  we  find 


In  "(r  I  V)    _  p 


n*(r  I  V)  =  J        ^"^^  I       -  "o^  ^ 


6w 


6n(r') 


or 

This  equation  we  will  call  the  Quantal  Hypernetted  Chain  or  QHNC 
equation. 

Following  similar  steps  with  the  functional 

u  _  n(r  I  U) 
^  "  n*(r  I  U) 

we  can  arrive  at 

or  what  shall  be  referred  to  as  the  QPY  equation. 

It  is  of  no  trivial  importance  that  the  QOZ,  QPY  and  QHNC 
equations  are  of  the  exact  same  form  as  their  classical  counterparts 
(this  in  addition  to  the  property  of  reducing  to  the  classical 
equation  in  the  limit).    This  is  because  phenomenologi cal  ways  of 
incorporating  the  effects  of  truncating  classical  functional 
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expansions  at  the  linear  term  may  still  be  of  use.    [In  the  classical 
HNC  case  such  a  truncation  is  equivalent  to  neglecting  bridge  diagrams 
in  the  Mayer  expansion  of  non-ideal  gases.]    This  will  be  considered 
further  in  Chapter  V. 

It  should  also  be  stressed  again  that  these  integral  equations, 
even  if  they  could  be  made  exact,  are  for  g^"(r)  =  n(r/v)/nQ  and 
not  the  pair  distribution  function.    This  is  evidenced  by  the  fact 
that  if  we  turn  off  interparticle  interactions,  V^"(r)  vanishes  and 
the  solutions  of  either  the  QHNC  or  the  QPY  equation  yield  g^"(r)  = 
1.    This  is  correct  for  the  density  distribution  but  the  fermion  pair 
distribution  has  an  exchange  hold.    This  difference  is  taken  up  in  the 
next  chapter. 

The  exact  correspondence  between  quantal  analogue  and  classical 
integral  equations  is  lost  when  one  considers  screening  effects  on 
long  range  potentials.    Rigorously  incorporating  screening  effects  is 
complicated  in  the  quantal  case  by  symmetry  effects.    It  can  be 
shown^^'^  that  the  best  single  particle  Hamiltonian  (best  in  the 
sense  that  the  statistical  average  over  the  square  of  the  difference 
of  the  exact  and  single  particle  Hamiltonian  is  minimized)  is  provided 
by  a  Hartree-Fock  potential,  not  the  Hartree.    On  the  other  hand  we 
will  see  in  the  next  section  that  as  far  as  density  distributions  are 
concerned  there  exists  an  exact  effective  potential  single  particle 
Hamiltonian.    In  the  quantal  integral  equation  approach  we  will 
consider  here  screening  of  the  reference  potential  will  ignore 
exchange  effects;  we  will  see  in  the  next  section  that  they  are 
actually  incorporated  in  the  closure  of  the  linear  functional 
expansion,  namely  the  QOZ  relation. 


-42- 


When  we  screen  the  reference  distribution  of  the  QPY  generating 
functional 


^  ^    n(?  I  U) 
n*(r  I  U^) 

with  the  Hartree  potential 


3U^(r)  =  (5U(r)  +  J  dr'(pU(r'  -  r)  (n(r'  I  U)  -  n^) 
its  linear  expansion  requires  the  variation 


6w 


6n(r') 


6(?  -  ?■)  _  1  '  V 

"o  "o    6n(r'  I  U) 


This  can  be  evaluated  by  invoking  the  chain  rule  over  the  Hartree 
potential,  which  in  Fourier  space  becomes 


6n*(r  I  U^) 
6n(r'  I  U) 


6n*(r  I  U^) 
6  -  pU^(r") 


S  -  3U^(r") 
■6n(r'  I  U) 


=  "oX*(0)  Fq{ 


6  -  3U^(r") 


6n(r'  I  U) 

Now  the  variation  of  the  Hartree  potential  yields 


6  -  3U^(r) 
6n(r') 


6  -  BU(r) 
6n(r') 


-  |iU(r  -  r') 


or 


6  -  3U^(r) 
Sn(r') 


}  = 


n^x(Q)  'a 


-  K{mr  -  ?')} 


from  which  a  few  simple  manipulations  yield 
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6w 


=  c^"(r-r')  +  Fq^  x*(0)  F{(iU(r-r')}  =  c^"(r-r')  +  M{|3U(r-r')} 
0  ^ 


5n(r') 

Plugging  into  the  functional  expansion  of  W[n]  we  obtain 


"^"^  '        =  1  +  ;  d?'(c^"(?  -  r')  (n(?'  I  v)  -  n  ) 
n*(?  I  V^)  0 

+  ;  dr'(M{flU(r  -  ?')})  (n(r'  I  v)  -  n^) 

In  terms  of  analogue  functions  and  employing  the  QOZ  relation  this 
becomes 

g    (r)  =  e  [1  +  Y    (r)  +  Fq  {x*(0)  F{n(?  I  v)  -  nj  F{(JV(?)}} 

=  e  [1  +  Y    (r)  +  M{;  dr'(pV(r  -  ?'))  (n(r'  I  V)  -  n^)}] 

=  e  [1  +  Y    (r)  +  M{3V^(r)  -  3V(?)}] 

This  equation  will  be  referred  to  as  the  QPYH  equation. 

Unlike  the  classical  case,  if  we  screen  the  QHNC  generating 
functional,  we  do  not  recover  the  QHNC  equation.  A  now  familiar 
procedure  yields  instead  the  result 

g^"(?)     e'  '^""^''^    ^^""^  ^  M{3V^(r)  -  pvc?)} 

to  be  known  as  the  QHNCH  equation.  We  see  that  for  both  the  QPYH  and 
QHNCH  equations  long  range  tails  cancel  but  the  potential  enters  in  a 
nontrivial  manner. 

Lastly,  in  the  classical  case  we  saw  how  the  HNC  equation  could 
also  be  derived  by  constructing  a  screened  potential  such  that  the 
first  variation  of  the  generating  functional  vanished.    The  same 
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procedure  can  be  applied  to  the  quantal  case.    Define  an  effective 
potential 

(5U.(r)  =  (5U(r)  -  J  dr'(c(r  -  r'))  (n(r'  I  U)  -  n^) 

I.  0 

(where  c(r)  is  the  quantal  pair  correlation  function  and  not  the 
analogue  quantity);  then  both  the  QHNC-  and  QPY-like  generating 
functional s 

W  =  In  '  U>  H^_ni^_LUi_ 

n*(r  I  U^)  n*(r  I  U^) 

have  vanishing  first  variations  and  both  yield  the  same  equation: 
an  - 

This  equation  is  distinct  from  the  QHNC  equation  but  like  the  QHNC 
equation  reduces  to  the  classical  HNC  equation  in  the  limit.    For  this 
reason  it  shall  be  termed  the  0HNC2  equation. 

Interestingly  enough  the  0HNC2  equation  can  be  obtained  from  an 
entirely  different  approach.    Instead  of  assuming  the  quantal 
bootstrap  relation  and  utilizing  Percus'  method  we  can  arrive  at  the 
QHNC2  result  directly  from  Kohn-Sham-Mermi n  density  functional 
theory. From  there  we  can  recover  the  quantal  bootstrap 
relation  (our  starting  ansatz)  and  so  justify  our  variant  equations  by 
making  what  we  will  see  as  a  reasonable  assumption  concerning  the 
effective  potentials  of  many  body  systems. 


•1 
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The  Density  Functional  Approach 

Following  the  zero  temperature  approach  of  Hohenberg  and 
Kohn^*^^  and  Kohn  and  Sham^^^  for  the  ground  state  properties  of  an 
electron  gas  in  an  external  potential  U(r),  Mermin  showed  that^^° 
i)  the  Helmholtz  free  energy  A  of  an  i nhomogeneous  electron  gas  in 
thermal  equilibrium  is  a  unique  functional  A[n]  of  the  one  electron 
density  n(r),  and  ii)  the  grand  potential 

Q[n]  =  A[n]  -  vi  J"  n(r)  dr 

is  minimum  for  the  equilibrium  density  n(r).    We  employ  the  grand 
canonical  rather  than  a  canonical  ensemble  because  macroscopic 
quantities  derived  from  a  non-interacting  Fermi  system  are  more  easily 
calculated  and  the  fundamental  variable,  namely  the  density 
distribution,  is  a  physical  observable  that  should  be  the  same  from 
either  ensemble.    The  free  energy  is  customarily  decomposed  as 

A[n]  =  1  J  v(r  -  ?')  n(r)  n(r')  dr  dr'  +  J  U(r)  n(r)  dr  +  A*[n]  +  A  [n] 

X  c 

to  separate  out  respectively  1)  solely  classical  electron-electron 
interaction  contributions,  denoted  by  v(r),  2)  the  external  potential 
contribution,  3)  the  contribution  A*[n]  that  a  non-interacting  Fermi 
system  with  density  n(r)  would  still  exhibit,  and  4)  all  other  effects 
are  grouped  to  define  the  unkown  exchange  and  correlation  functional 
A^^[n].    The  minimum  property  of  the  grand  potential  leads  to  the 
functional  equation 

6A*[n] 

— +  U  f^(r)  -  }i  =  0  (3.5) 
6n(r) 
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where  u  is  the  chemical  potential  and  U^^^d)  is  defined  by 

SA  rn] 

Ug^^(r)  =  U(r)  +  J  dr  V(r  -  ?')  n(r)  + 

6n(r) 

Equation  3.5  is  equivalent  to  the  equation  for  non-interacting 
electrons  in  this  effective  potential;  that  is  we  have  n(rlU)  = 
n*(rlUg^^).    The  non-interacting  fermion  distribution  n*(rlU)  is 
given  by  the  following  system: 

(-  ^  A  +  Ug^^)       =  e.  fi  =  M  =  1  units 

n*(r  I  Ug^^)  =  I  f(e.)  I  ^.(r)  1^ 

where  the  sum  over  the  index  "i"  refers  to  all  bound  and  continuum 
eigenstates  and  f(e)  is  the  Fermi-statistics  occupation  number 

1-1 


f( 


e)  =  {1  +  exp[p(e.  -  li)]}' 


dependent  on  the  temperature  and  the  chemical  potential.    The  latter 
is  determined  as  usual  from  the  total  number  of  electrons  by  spatially 
integrating  over  n*(r). 

For  the  jellium  model  of  the  electron  gas  the  external  potential 
must  not  only  include  the  source  of  the  external  potential  but 
the  effect  of  the  neutralizing  background  as  well 

U(r)  =  U^Cr)  -  J  dr'  n^  V(r  -  ?') 
so  that  the  effective  potential  is 

Ueff(r)  =  U^Cr)  +  ;  dr'  V(r  -  r')  (n(?  I  U^)  -  n^  + 
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At  large  radii,  the  displaced  electron  charge  vanishes,  and  any  long 
range  tail  of  the  source  potential  is  cancelled  on  the  assumption  of 
perfect  screening,  so  that  the  assymptotic  value  of  the  effective 
potential  is 


6A 


W^^">=  6n 


xc 


Applying  eq.  5  at  large  radius  gives 
6A*  6A. 


6n. 


=  ^  - 


xc 


6n 


or,  equivalently 


v2  „-  3/2 


n„  =  —  3    ^  I^/^Pii) 


0  2 
ir 


where  the  standard  definition  of  the  Fermi-Dirac  functions  has  been 
used: 


I  (n)  E  J  y* 

«  0  1.  e^-^ 


It  is  convenient  to  shift  the  zero  point  of  our  energy  scale  to  define 
a  new  effective  potential  which  goes  to  zero  at  infinity 

Ugff(r)  =  Uq(?)  +  J  d?'  V(?  -  r)  (n(r  I  U  )  -  n  ) 


6A 


+  (- 


xc 


6n 


6A. 


)  -  (- 


xc 


6n 


) 


Up  to  this  point  our  equations  are  exact;  the  exchange 
correlation  part  of  the  free  energy  functional  is,  however,  quite 
unknown.    But  if  we  make  a  functional  expansion  to  first  order 
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6A. 


xc 


6n 


6A 


)  =  (- 


xc 


5n 


)  +  ;  K(r  -  r')  (n(r'  I  U^)  -  n^)  dr' 


5  A 


K(r  -  r')  E  (- 


xc 


6n(r)  5n(r') 


we  obtain  an  approximate  integral  equation  in  terms  of  the  unknown 
kernal  K(lr-r'l): 


n(?  I  V^)  -  n^  =  n*(-r  I  U^^^)  -  n^ 


^eff  =  ^0     -f"  dr'{V(r  -  r')  +  K(r  -  ?■)}  {n(?'  I  U^)  -  n}^ 

Since  the  unknown  kernal  is  independent  of  the  external  potential 
U(r),  we  can  determine  this  function  by  considering  the  case  of  a  very 
weak  external  potential  where  the  linear  response  formula  can  be 
employed  on  both  sides  of  the  first  of  the  above  equations  with 
impunity.    The  left  hand  side  yields 

6n(0)  E  FQ{n(?  I  U^)  -  n^}  ^  -  n^p  x(Q)  U^CQ) 
while  the  right  hand  side  gives 

FQ{n*(r  I  Ug^^)  -  n^}  :  -  n^3x*(Q)  U^^^CQ) 

=  -  n^^x*(Q)  {Uq(0)  +  [V(0)  +  K(0)]  8n(Q)} 

Using  the  last  equation  we  can  show  that  the  Fourier  transform  of  the 
unknown  function  K(r)  is 


K(0)  =  (- 


V    x(0)  x*(0) 


)  -  V(Q) 
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Due  to  this  assignment  the  effective  potential  reduces  in  form  to  the 
previously  defined  screened  potential 

Uj.(?)  =  Uq(?)  -  ;  d?'  (1  c(?  -?')  (n(r  I  U^)  -  n^) 
and  we  obtain  the  relation 
n(r  I  U^)  =  n*(r  I  U^) 

Note  that  the  source  external  potential  is  completely  arbitrary,  and 
so  this  result  is  a  generalization  of  the  QHNC2  equation,  which 
follows  immediately  when  the  external  potential  is  taken  to  be  the 
interparticle  interaction.    From  this  generalized  result  we  can  work 
backwards  and  obtain  the  quantal  bootstrap  relation,  and  so  justify 
the  other  (QPY,  QHNC,  etc)  equations  we  have  obtained. 

This  is  done  by  now  making  the  assumption  that  the  average  one 
body  potential  felt  by  a  constituent  electron  due  to  a  fixed  test 
charge  6e,  namely 

Vt  =Vtest   -  ^  d-r'(^  c(-r  -  r'))  (nCr'  I  V^^^^)  -  n^) 

where 


is  equal  to  the  average  one  body  potential  felt  by  a  free  test  charge  Se, 

arising  from  a  fixed  electron  in  an  electron  gas,  namely 

V^-^g  =  (5e)  l^(r) 

W(r)  =  4Tr{  e8(r)  +  e[n(r  I  U  =  e^/r  -  n  ]} 
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In  the  limit  of  small  Se  we  can  again  invoke  linear  response  with 
impunity  and  solve  for 

F^{n(r  I  U  =  V)  -  n^}  =  -  1 

thus  recovering  the  quanta!  bootstrap  relation. 

He  see  that  our  density  functional  approach,  aside  from 
justifying  the  quanta!  bootstrap  ansatz  of  the  Percus  method,  gives 
physical  insight  to  its  meaning.    Even  though  it  is  not  an  exact 
approach,  it  clearly  points  out  its  limitation,  as  we  approximated  the 
variation  of  the  exchange-correlation  free  energy  functional  by  a 
first  order  functional  Taylor  expansion.    This  type  of  truncation  is  a 
property  shared  by  the  Percus  method  and  ultimately  may  be  responsible 
for  the  self-consistency  of  the  two  approaches. 


CHAPTER  IV 
THE  BASIC  QHNC  EQUATION 


An  alternative  ansatz  for  the  extension  of  classical  integral 
equations  into  the  quantal  regime  lies  in  the  construction  of 
effective  pair  potentials.    These  incorporate  the  bulk  of  quantal 
effects;  for  example  in  mul ti component  plasmas  quantal  effects  prevent 
the  collapse  of  particles  in  an  attractive  coulomb  field.  More 
generally  the  effective  potential  yields  finite  values  for  the  pair 
distribution  factor  g(r)  in  the  limit  of  zero  interparticle  separation 
r. 

Such  classical  pseudopotential s  have  been  constructed  by 
requesting  an  approximate  account  of  the  N-body  Slater  sum^^^~^^^ 
or,  in  an  approach  which  can  in  principle  be  implemented  exactly  but 
is  more  limited  in  its  physical  range  of  applicability,  by  reproducing 
the  two-body  Slater  sum^^^"^^^ 

-  (r)  -  (JH^^^ 

e  =  (A-^'^^/my'^  <r  I  e      ^     I  r>  (4.1) 

"2'^  =  -  ^  v2  +  V(r) 

It  can  be  shown  that  this  Binary  Slater  Sum  approximation  correctly 
represents  the  derivative  of  the  classical  pseudopotential  for  the 
N-body  problem  in  the  limit  r  tends  to  zero.^^^ 
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As  shown  in  the  previous  chapter,  the  QHNC  equation  for  n(rlv)  is 
also  equivalent  to  the  classical  HNC  equation  with  an  effective 
potential  defined  by 

-  P^eff^^^     n*(r  I  V) 
e  =  — — — 

"o 

where  n*(r|v)  is  the  density  distribution  of  non-interacting  fermions 

in  the  presence  of  an  external  potential  of  the  form  of  the 

interparticle  interaction—the  bare  coulomb  potential.    Normally  this 

entails  the  solution  of  the  Kohn-Sham-Mermin  equations That 

is  calculating  n*(r|v)  is  the  solution  of  the  following 
system:120-121 

(-  ^  v2  .  V)  ^.  =  c.  f. 

n*(r)  =  I  f(e.)  (4.2) 
i  111 

where  f(e.)  is  the  Fermi-statistics  occupation  number 
f(e.)  =  {1  +  exp[3(c.  -  \i)]y^ 

depending  on  the  temperature  and  the  chemical  potential  fixed  by  the 
total  number  of  particles.    The  index  "i"  runs  over  all  continuum  and 
bound  (if  existent)  states. 

The  treatment  of  the  delocalized  states  proposed  by  Friedel  for 
impurities  at  zero  temperature      can  be  applied  to  the  nonzero 
temperature  case  too.    The  eigenfunctions  ■i^^^^  of  eq.  4.2  for 
energy 
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are 

(|>j^^(r)  is  normalized  by  requiring  that  it  have  the  asymptotic  form 
^^^■^(.r)  -  ^  sin[Kr  -  e  J  +  n^CK)] 

n-jCK)  being  the  phase  shifts  in  the  potential  V(r)  and  n^(K  =  «>)  =  0. 
Overall  normalization  is  accomplished  by  requiring  that  the 
wavefunctions  vanish  at  a  spherical  wall  of  very  large  radius  R.  We 
assume  that  R  is  so  large  that  most  of  the  integrand  is  approximated 
by  its  asymptotic  form: 

1  =  I  A^^ 1^  ;   dr  sin^CKr  -  e  J  +  n^] 

with  the  boundary  conditions 

KR  +  T]^iK)  =  e  J  =  mir  m  =  0,  1,  2,  3.  .... 

incrementing  m  by  unity  for  a  given  L  yields  the  density  of  states 
(including  spin  multiplicity)  as 

9i(K)  =f  ^1  ^  rIkHiCK)) 

while  the  normalization  factor  for  very  large  R  is  given  by 
A^^  =  v27r  (1  -0(^)) 

(These  same  results  can  be  derived  for  a  long  range  coulombic 
potential .^^^  36.23^    ^  straightforward  derivation  gives  for 

the  density 

n*(r)  -       =  I  f<ic^)  %  %  ^       !    f(^)  dK  I  (2e  +  1)  [<|,2    -  2^(Kr)l 
b  ir     0  1  Ml 

! 

I 

■I 
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where  subtracting  off  the  uniform  density  result  accelerates 
convergence  of  the  integral  and  summation.    We  note  that  in  general  a 
potential  may  support  a  number  of  bound  states  and  we  have  included 
them  al so. 

The  OHNC  Effective  Potential 

We  thus  see  that  the  QHNC  pseudopotential  is  essentially  a 
one-particle  Schrodinger  equation  calculation.    This  is  a  feature  in 
common  with  the  Binary  Slater  Sum.    It  is  important  to  contrast  the 
two  however.    The  Binary  Slater  sum  (BSS)  is  used  to  calculate  the 
pair  distribution  function  g(r)  while  the  QHNC  equation  yields 
n*(r|v);  this  difference  is  important  in  consideration  of  quantum 
effects.    It  should  be  noted  that  in  the  BSS  approach  the  Schrodinger 
equation  is  solved  with  the  use  of  the  reduced  mass.    The  matrix 
element  in  eq.  4.1  can  be  evaluated  by  inserting  a  phvsi cal  basis  set 
of  properly  symmetrized  or  anti symmetrized  wavefunctions ,  depending  on  i 
the  statistics  of  the  two  particles.    Parallel  spin  versus 
antiparallel  spin  contributions  can  be  calculated  separately.    In  the 
BSS  approach  the  effective  potentials  are  density  independent.    In  the 
QHNC  approach  the  Kohn-Sham  system  of  equations  can  also  be  written  in 
the  form  of  a  matrix  element 

n*(r  I  V)  =  <r  I  f  I  r> 
where  the  operator  "f"  is  defined  as 

f  =  {n  +  exp  m  -  p]}"^  - 
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H  is  the  coulomb  hamiltonian  operator,  mu  the  chemical  potential  and 
we  have  included  a  variable  "eta"  to  be  set  equal  to  unity  to  describe 
Fermi  statistics.    In  the  QHNC  approach  the  external  potential  has  no 
spin  attribute,  and  the  matrix  element  is  evaluated  using  a  complete 
basis  set  with  no  regard  to  statistics;  the  statistics  are  built  into 
the  functional  form  of  the  matrix  element.    The  QHNC  effective 
potentials  are  implicitly  dependent  on  the  density  through  the 
chemical  potential. 

All  these  subtle  differences  are  manifestations  of  the  fact  that 
in  the  QHNC  approach  g(r)  must  be  obtained  from  n*(r|v)  by  a  further 
calculation.    One  obvious  advantage  into  this  splitting  of  tasks  is 
that  the  BSS  pseudopotential  breaks  down  at  zero  temperature,  whereas 
the  QHNC  remains  valid. 

A  surprising  feature  of  the  QHNC  effective  potential  is  that  it 
is  analytically  calculable.    This  allows  one  to  bypass  the 
computationally  time  consuming  Kohn-Sham-Mermin  system  of  equations. 
The  starting  point  for  the  solution  is  the  fact  that  the  effective 
potential  is  related  to  the  r^  =      limit  of  the  off  diagonal 
density  matrix  for  fermions  in  an  external  coulomb  potential : 

p(f^  I  r^)  =  <f^  I  f  I  r^> 

Spherical  symmetry  reduces  the  functional  dependance  of  the  off 
diagonal  density  matrix  does  from  the  full  six  variables  of  r^  and 
12  to  only  the  magnitudes  r^ ,      and  the  angle  between  r^  and 
T^.    This  is  easily  seen  if  we  expand  the  matrix  element  using  a 
complete  set  of  spherical  coulomb  wavefunctions : 
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dL  c     ri^2     2m          ,  (1  +  1),  ^  « 
X  S  +  [K    -  ~^  ~  -  I   0^  S  =  0 

Using  the  addition  theorem  of  spherical  harmonics  then  yields 
where  the  1-wave  radial  density  matrix  components  are  defined  by 

00 

Pl(r^  I  r2)  =  ;    (spin  J  dK)  — ^-^          S*^(r^)  S^^^r^) 

^  ^  2m  ^  ^ 

Furthermore,  because  the  potential  is  coulombic,  there  is  an 

additional  symmetry:    the  matrix  element  "f"  also  commutes  with  the 

1 24 

Runge-Lenz  vector.        This  reduces  the  functional  dependance  of  the 
off  diagonal  density  matrix  to  the  two  variables  "x"  and  "y": 

7         7  ^  17 

X  =  (r^  +  r2)/2  y  =  (r^  +  r2  -  2r^  cose)"^/2 

The  solution  then  follows  the  steps  of  Hostler^^^"^^^  and  Storer^^^ 
except  that  the  Bloch  equation  satisfied  by  the  off-diagonal  canonical 
density  matrix  is  replaced  by  a  Fermi  statistics  generalization 

[(H  -  1^)  n  +  n       +  fj^]  p(r^  I  r2)  =  0 

(Here  H  acts  only  on  the  r^  coordinates  and  we  set  eta  equal  to  unity 
at  the  end  of  the  calculation.)    Because  of  the  aforementioned 
symmetries,  the  solution  to  this  equation  reduces  to  a  functional  form 
identical  with  the  s-wave  radial  density  matrix  component  (however  in 
the  variables  x  and  y): 
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The  diagonal  limit  r]  =  r2  =  r  is  given  by 
n*(r)  .  p(?  1  7)  =  -  ^         p  (X  +  y)  x 

°^  ay 


-  y) 


y=o 


The  second  derivatives  of  the  radial  coulomb  wavefunctions  can  be 
eliminated  in  favor  of  the  radial  wavefunctions  themselves  by 
employing  the  radial  coulomb  differential  equation.    We  thus  obtain 


n*(r)  =  ^  J  dK  ( 


1 


2ir' 


J(K,  r)  =  { 


3  rf<V  ,  , 
e    t  0^   -  >i]  +  1 


)  J(K.  r) 


d?  ^Ko^^^ 


2m 
2 


-  -  K^)  I  S^^(r)|2} 


Numerical  Calculations 

The  effective  potential,  obtained  from  the  ideal  fermion  density 
distribution  in  a  coulomb  potential,  will  be  calculated  from  the 
analytic  formula  derived  above.    For  comparison  with  the  classical  HNC 
equation,  we  will  scale  lengths  in  units  of  the  interparticle 
separation  x  =  r/a  where 

n^  =      u  a3)-l 

Temperature/density  points  will  be  expressed  by  the  classical  coupling 
parameter 
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along  with  a  quanticity  temperature  (two  divided  by  the  temperature  in 
rydbergs) 

13*  E  fime^/fi^ 

In  the  classical  limit  all  functions  are  solely  dependent  on  the 
parameter  r^;  we  expect  any  dependance  on  the  quanticity 
temperature  to  indicate  the  pervasiveness  of  quantum  effects. 

The  calculations  were  performed  corresponding  to  the  temperature 
density  points  of  tabulated  results  provided  by  Pokrant     for  later 
comparison.    These  consisted  of  values  at  constant  density,  as 
measured  by  the  dimensionl ess  ground  state  coupling  parameter 

r^  =  a/a^ 
s  0 

where  a^  is  the  bohr  radius,  as  well  as  a  dimensionl ess  degeneracy 
temperature 

t  =  Kt/Ep 

where      is  the  ground  state  Fermi  energy.    The  results  are 
presented  in  the  following  graphs. 

In  Fig.  4.1  we  present  plots  of  the  ideal  fermion  density  in  the 
presence  of  an  external  coulomb  potential.    For  comparison  these  are 
plotted  together  with  two  approximations.    The  first  is  merely  the 
classical  limit 

n*(r  I  V) 

The  second  is  the  Thomas-Fermi  approximation. It  is  obtained  by 
expanding  the  diagonal  density  matrix  element  "f"  using  a  complete  set 
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of  plane  waves.    If  we  then  neglect  the  non-commutativity  of  the 

potential  and  kinetic  energy  operators  in  the  exponential  (the 

neglected  terms  can  be  systematically  calculated  by  Trotter 
1 23 

formulas      and  shown  to  vanish  for  large  r)  we  then  obtain 

'j 

n*(r  I  V)     ^1/2  (Pv^-  r/x)  I 
"o        "  ^1/2^^^^ 

where  the  Fermi-Dirac  integral  of  order  p  is  defined  by  the  integral 

0      1  +  e-' 

(In  this  formula  the  symbol  r  denotes  the  gamma  function.)    It  is 
seen  that  the  exact  result  quickly  merges  with  the  Thomas-Fermi 
approximation. 

By  using  the  relation 

dF  (X) 

we  find  that  asymptotically  the  Thomas-Fermi  approximation  has  the  form 

F   ^  (3n) 
n*(r  I  V)  ~  1     I     "  2 

^        1  ^^V"^ 
^  2 

This  is  indeed  the  correct  asymptote  of  the  analytic  solution  as  can 
be  verified  by  looking  at  the  appropriate  expansion  of  the  spherical 

1 29 

coulomb  wavefunctions.        This  tells  us  that  the  effective 
potential  in  the  QHNC  approximation  has  a  long  range  1/r  tail—but 
with  a  strength  modified  from  the  classical  theory.    The  QHNC  tail 
contains  the  additional  factor 
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F    ^  ((5p)/F  ^ 
"2  "^2 

not  found  in  the  Binary  Slater  Sum  approach.    In  Fig.  4.2  we  present 
plots  for  the  QHNC  effective  potential.    For  comparison  the 
corresponding  effective  potentials  of  the  Binary  Slater  Sum 
approximation  are  presented  in  Fig.  4.3. 

The  approach  to  solving  the  HNC  system  of  equations  with  long 
tail  potentials  is  well  known — we  follow  the  convention  of  Ng^^^  and 
renormalize  the  constituent  functions  by  subtracting  out  an 
analytically  transformable  long  range  function  of  the  form 

^  erFc(ar) 

We  employed  the  Ng  method  of  solution  (essentially  a  picard  iteration 

with  an  accelerated  convergence  procedure)  although  it  should  be 

mentioned  that  two  new  algorithms  have  recently  been  published^^^"^^^ 

that  employ  a  hybrid  Picard/Newton-Raphson  iteration  procedure  which 

promises  faster  convergence. 

The  results  of  the  integral  equation  solver,  that  is  the  analogue 

pair  distribution  function  n(.r/yj)/r\^,  are  presented  in  Figs.  4.4-4.8. 

This  distribution  is  of  theoretical  interest  in  its  own  right,  for 

example  in  the  study  of  impurities  in  metals  or  positron  annihilation 
1 32 

rates.        However  we  will  limit  our  concern  to  the  approximate 
calculation  of  the  pair  distribution  function.    A  method  of  calculating 
g(r)  based  in  part  on  the  information  contained  in  n(r|v)/n  ,  and  a 

0 

comparison  of  results  with  the  Binary  Slater  Sum  approximation  and 
that  of  Pokrant  will  be  taken  up  in  the  next  two  chapters. 
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Figure  4.1 

Plots  of  n*(r|v)/no-l  versus  x  =  r/a  for  a  density  of  rg  =  2 
at  temperatures  of  kt/Ef  of  .1,  .5,  1.0,  2.0,  20.0.  Storer  denotes 
the  analytic  solution,  following  Russian  literature  the  Thomas-Fermi 
approximation  has  been  denoted  Hartree,  and  Boltzmann  is  the  classical 
approximation. 


tau-  8.10 


0.  -f- 


T 


1.0 


 1 — 

,5  2.0 
x-r/a 


2.5 


3.1 


3.5 


Figure  4.2 

Plots  of  the  QHNC  effective  potential  versus  x  =  r/a  for  a 
density  of  rg  »  2  at  temperatures  of  kt/Ef  of  .1,  .5,  1.,  2.,  20. 
We  have  plotted  the  potential  in  the  form 

I  Veff('^> 

A  purely  classical  coulombic  interaction  would  correspond  to  a  constant 
values  of  unity.    Storer  represents  the  analytic  solution;  Hartree  is 
the  Russian  literature  denotation  of  the  Thomas-Fermi  Approximation, 
which  is  truncated  near  the  origin  to  prevent  divergences. 
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Figure  4,3 

Plots  of  the  Binary  Slater  Sum  effective  potential  versus  x  = 
for  a  density  of  rj  =  2  at  temperatures  of  kt/Ef  of  .1,  .5,  1., 
2.,  3.,  5.    We  have  plotted  the  potential  in  the  form 

I  Veff<K) 

A  purely  classical  coulombic  interaction  would  correspond  to  a 
constant  values  of  unity.    Note  that  as  we  approach  classical 
temperatures  the  curve  reaches  this  constant  value  quickly,  and  in 
fact  first  overshoots  this  value. 
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Figure  4.4 

Plots  of  density  times  the  Fourier  transform  of  h(r)  =  g(r)  -1 
versus  k  =  q*a  (q  the  reciprocal  space  to  r)  for  a  density  of  rc  -  2 
at  temperatures  of  kt/Ef  of  .1.    Free  denotes  the  ideal  fermion  pair 
distribution  function,  HNC  the  results  of  the  QHNC  calculation  where 
g(r)  represents  the  analogue  pair  distribution  n(r/v)/no.  Chihara 
denotes  the  results  of  a  yet  to  be  discussed  construction  of  g(r)  from 
the  analogue  g(r). 
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Figure  4.5 

Plots  of  density  times  the  Fourier  transform  of  h(r)  =  g(r)  -1 
versus  k  =  q*a  (q  the  reciprocal  space  to  r)  for  a  density  of  rj  =  2 
at  temperatures  of  kt/Ef  of  0.5.    Free  denotes  the  ideal  fermion 
pair  distribution  function,  HNC  the  results  of  the  QHNC  calculation 
where  g(r)  represents  the  analogue  pair  distribution  n(r/v)/no. 
Chihara  denotes  the  results  of  a  yet  to  be  discussed  construction  of 
g(r)  from  the  analogue  g(r). 
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Figure  4.6 

Plots  of  density  times  the  Fourier  transform  of  h(r)  =  g(r)  -1 
versus  k  =  q*a  (q  the  reciprocal  space  to  r)  for  a  density  of  rg  =  2 
at  temperatures  of  kt/Ef  of  1.0.    Free  denotes  the  ideal  fermion 
pair  distribution  function,  HNC  the  results  of  the  QHNC  calculation 
where  g(r)  represents  the  analogue  pair  distribution  n(r/v)/no. 
Chihara  denotes  the  results  of  a  yet  to  be  discussed  construction  of 
g(r)  from  the  analogue  g(r). 
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Figure  4.7 

Plots  of  density  times  the  Fourier  transform  of  h(r)  =  g(r)  -1 
versus  k  =  q*a  (q  the  reciprocal  space  to  r)  for  a  density  of  rg  »  2 
at  temperatures  of  kt/Ef  of  2.0.    Free  denotes  the  ideal  fermion 
pair  distribution  function,  HNC  the  results  of  the  QHNC  calculation 
where  g(r)  represents  the  analogue  pair  distribution  n(r/v)/no. 
Chihara  denotes  the  results  of  a  yet  to  be  discussed  construction  of 
g(r)  from  the  analogue  g(r). 
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Figure  4.8 

Plots  of  density  times  the  Fourier  transform  of  h(r)  =  g(r)  -1 
versus  k  =  q*a  (q  the  reciprocal  space  to  r)  for  a  density  of  rj  =  2 
at  temperatures  of  kt/Ef  of  8.0.    Free  denotes  the  ideal  fermion 
pair  distribution  function,  HNC  the  results  of  the  QHNC  calculation 
where  g(r)  represents  the  analogue  pair  distribution  n(r/v)/no. 
Chihara  denotes  the  results  of  a  yet  to  be  discussed  construction  of 
g(r)  from  the  analogue  g(r). 


CHAPTER  V 
LOCAL  FIELD  CORRECTIONS 


We  have  seen  in  the  previous  chapter  how  classical  integral 
equations  generalized  into  the  quantal  region  can  yield  accurate 
results  for  the  quantity  n(rlv),  namely  the  density  distribution  of 
interacting  particles  in  the  presence  of  an  external  potential  of  the 
form  of  the  interparticle  interaction.    This  in  effect  is  equivalent 
to  knowledge  of  the  static  response  function  x(k).    If  from  this 
information  the  dynamic  response  x(k.,w)  could  be  extracted,  then  we 
could  easily  calculate  the  static  structure  factor  or  pair 
distribution  function  via  the  fluctuation  dissipation  theorem 

S(K)  =  ;     S(Kco)  d<o  =  i  ;     -rrr          Im  (Jx(Kco)  dco 

These  operations  would  constitute  the  quantal  generalization  of  the 
classical  bootstrap  assignment  g(r)  =  n(r|v)/nQ. 

To  demonstrate  that  such  a  program  is  plausible,  we  note  that  the 
analytical  properties  of  the  dynamic  susceptibility  make  it  possible 
to  write  an  exact  "dispersion"  relation  for  x"^(k,z)  of  the  form^^"^ 

X'^KZ)  =  x"^K)  +  ^  [Z^  +  iZK^  0(KZ)] 
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This  particular  relation,  which  is  one  of  many  possible 
representations,  introduces  an  unknown  function  e(k,z)  which  is 
analytic  in  the  upper  half  plane.    The  physically  meaningful  dynamic 
response  is  given  in  the  limit  z  =  w  +  ie,  and  in  this  limit  the  real 
and  imaginary  parts  of 

G(K,  Z  =  to  +  ie)  =  0'(Kco)  +  ie"(K(o) 

are,  like  the  dynamic  susceptibility  itself,  related  by  the 
Kramers-Kronig  relations.    The  function  e'(k,w)  is  known  to  be  a 
positive,  even  function  of  frequency,  as  follows  from  the  known  parity 
of  x'(k,w)  and  x"(k,w)  and  the  fact  that  S(k,w)>  0.    The  utility  of  the 
function  e(k,z)  is  that  it  interpolates  smoothly  between  the  known 
results  for  x'(k,z)  in  the  limits  of  large  and  small  z,  and  being 


less  sensitive  than  x(k,z)  itself,  is  more  amenable  to  approximation. 

By  subtracting  out  a  corresponding  expression  for  a  noninteracting 
Fermi  system  we  can  express  the  dynamic  susceptibility  as 


Here  we  have  grouped  all  of  the  actual  physics  into  a  dynamic 
effective  potential  Vg^^(k,z),  actually  a  functional  of  the  static 
response  x(k)  and  the  unknown  function  D(k,z).    We  will  show  through 
the  use  of  the  Mori  projection  operator  technique^^^  that  under 
reasonable  assumptions  the  dynamic  effective  potential  can  be 
approximated  by  the  static  (z  =  0)  limit  of  the  above  equation 


x(KZ)  = 


(5.1) 


1  -  V 


eff 


(KZ)  x  (KZ) 


V 


eff 


(KZ)  =  -  ( 


X(K)  Xq(K) 


) 
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It  should  be  noted  that  although  x^{k,z)  in  eq.  5.1  is 

usually  taken  to  be  the  free  response  function,  this  choice  is  not  a 

unique  one;  x^Ck.z)  plays  a  role  of  a  reference  system  in  terms  of 

which  the  expression  for  x(k,z)  is  constructed.  A  possible 

improvement  in  the  procedure  may  be  achieved  by  using  an  adjusted 

choice  of  y.^(k,z)  that  includes  self-energy  effects,  for 
14-19  137 

example.        '        This  method  of  improvement  will  not  be  taken  up 
in  this  thesis  in  lieu  of  other  avenues  of  extension.    We  will  retain 
the  free  response  as  our  reference  system;  this  will  allow  us 
comparison  with  certain  mean  field  theories. 


Mean  Field  Theories 


Mean  Field  Theories  bypass  the  dependence  of  v^^^Ck.z)  on 
x(k)  and  directly  try  to  relate  V^^^  to  the  bare  interparticle 
potential  by  solving  approximately  the  many  body  problem.  (Something 
we  have  already  done  via  our  integral  equations.)    The  efective 
potential  is  usually  taken  to  be  "static,"  that  is  independent  of  z, 
and  is  commonly  expressed  in  the  form 

Vg^^(K)  =  Vq(K)  [1  -  G(K)] 

where  V^(k)  is  the  Fourier  transform  of  the  bare  interparticle 
potential  and  G(k)  is  called  the  local  field  correction  factor.  A 
review  of  such  local  field  function  theories  can  be  found  in 

1  og 

Kugler.        The  many  body  part  of  the  problem  is  generally  solved  by 
the  Vlasov  kinetic  equation      or,  what  amounts  to  the  same  thing, 
appropriately  decoupling  the  equations  of  motion  for  the  density 


-72- 


fluctuations J'*^   Whatever  route,  at  some  point  collisions  between 

particles  are  neglected  and  the  particles  are  assumed  to  move  in  a 

mean  or  average  field  based  on  the  collective  motions  of  all  the  other 

particles.    The  physical  significance  of  the  correction  factor  is  thus 

to  incorporate  some  of  the  effects  of  short  range  correlations  between 

individual  particles  not  taken  care  of  by  the  average  field. 

Mean  field  theories  for  the  most  part  have  been  limited  to  the 

ground  state  quanta!  electron  gas.    These  have  disagreed  on  the  form 

of  the  local  field  correction  factor  on  three  main  points. 

141_145 

First  some  workers  have  been  led  to  the  conclusion  that 

G(k)  is  a  universal  function  of  k/k^  (k^  the  Fermi  wave  vector), 
146_i4g  29-31 

while  others  '         contend  that  G(k)  is  density  dependent. 

Second,  there  is  wide  diversity  of  opinion  regarding  the  value  of 

1 50-1 51 

G(k)  in  the  limit  of  large  wavevector  k.  At  this  limit,  a 

constant  1/3,  irrespective  of  density,  has  been  given  by  Geldart  and 
Taylor. Rajagopal^^^  and  others, ^^^'^^^"^^^  while  that  of 

Togio  and  Woodruff^^^  and  Kugler^^^  use  2/3.    On  the  other  hand 

29  30 

Singwi  et  al .    '     derive  this  value  to  be  l-g(O)  in  terms  of  the 

pair  distribution  function  g(r)  at  the  origin,  and  that  of  Vashista 
31 

and  Singwi  is 

1  -  g(o)  -  ap(^) 

(where  a  is  an  adjustable  parameter  determined  self-consistently)  in 
contrast  with  the  result  of  Niklasson^'*^  who  gives  2/3(1  -  g(0))  to 
this  limit. 

Thirdly,  there  are  several  authors^^^'^^^"^^^'^^^  who  have 
shown  the  local  field  factor  G(k)  has  a  maximum  around  k  =  2k^. 


-73- 


Kugler      among  them  further  indicates  the  slope  of  G(k)  has  a 

logarithmic  singularity  at  k  =  2k^  in  addition  to  a  maximum.    On  the 

29-31 

contrary,  G(k)  of  Singwi  and  collaborators         has  neither  maximum 

nor  singularity  in  its  slope  around  k  =  2k^.    It  should  be  noted 

that  the  behavior  of  G(k)  in  the  large  wave  vector  region  is  quite 

extremely  sensitive  to  small  differences  in  x(k)  so  that  it  offers  a 

stringent  test  for  the  distinction  of  various  theories. 

In  the  case  of  the  degenerate  electron  gas,  there  is  neither 

computer  simulation  nor  experiment  to  be  used  as  a  criterion  to  check 

these  diverse  conclusions  concerning  G(k).    However,  at  present  the 

29-30 

prescription  of  Singwi  et  al .         is  regarded  as  fairly  successful 

in  treating  a  degenerate  electron  gas  and  is  relatively  easy  to 

program  numerically.    It  encompasses  in  the  classical  limit,  as  do  the 

integral  equation  methods  of  Chapter  II,  the  Hyper-Netted  Chain  ' 

equation,  which  is  known  to  give  the  best  fit  to  the  pair  distribution 

function  of  computer  simulations.    For  this  reason  the  derivation  of 

the  Local  field  correction  factor  in  the  STLS  theory  will  be  presented 

here  as  a  basis  for  later  comparisons. 

The  STLS  Method 

o 

Following  Singwi  we  will  derive  the  density  response  function 
by  following  the  time  evolution  of  the  off  diagonal  single  particle 
density  function 

P  (X  I  X' ;  t)  E  <  ^■^(xt)  t  (x't)> 
a  a  CT 

i 
1 
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If  we  employ  a  Heisenberg  equation  of  motion  approach  with  a 
HamiUonian  containing  an  external  potential  Vg^^Cx.t),  we  obtain 

3t  "  2m  ^-  "  ^H^^^^  "  V^Cx't)}  p^(x  I  x';  t)  = 

-  r  d     [V(x  -  X")  -  V(x'  -  X")]  <^''^(xt)  n(x"t)  t  (x'  t)> 

where  the  Hartree  potential 

V^(xt)  E  Vg^^(xt)  +  ;  dx'  V^(x  -  x")  p(x".  t) 

depends  on  the  i nterparti cl e  interaction  potential       and  the 
average  single  particle  density  distribution 

p(xt)  =  <n(xt)>  E  <  T  'f'^Cxt)  t  (xt)> 

S 

We  were  able  to  regroup  the  terms  of  the  equation  of  motion  into  its 
present  form  by  defining  the  cumulant  bracket 

<'f'^(xt)  n(x"t)  ^  (x't)>  E 
a  a  C 

<?'^(xt)  n(x"t)  t  (xt)>  -  <n(x"t)>  <^'"^(xt)  (x't)> 
a  a  a  a 

Our  equation  of  motion  is  equivalent  to  an  infinite  set  of  equations 
of  motion  for  observable  physical  quantities  with  classical 
equivalents.    This  follows  from  the  fact  that  the  Wigner  distribution 
function^^'^^^ 

fp^(Rt)  =  ;  d?  e"  '^^'^'^  <<tl  (R  -  ^  ?.  t)  t^(R  +  1  ?,  t)> 
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has  properties  analogous  to  the  classical  phase  space  distribution 
f(r,£;t).    In  particular  it  can  be  used  to  construct  physical 
observables  such  as  the  single  particle  density  distribution 

p(Rt)  =  I    f  (Rt) 
po 

the  particle  current  density 

J(Rt)  =  m"^  ^    P  f  (Rt) 
^  pa 
pa 

and  the  kinetic  stress  tensor 

r(Rt)  =  m"M    P  P  f  (Rt) 

pa 

pa  ^ 

By  expanding  the  equation  of  motion  for  the  off  diagonal  density 
distribution  about  its  diagonal,  that  is  in  powers  of  r  =  x'  -  x,  the 
coefficients  of  the  first  few  powers  of  r  yield  the  beginning  of  a 
hierarchy  of  one  particle  equations.    First  we  have  the  usual 
continuity  equation 

1^  p(Rt)  +  V«J  (Rt)  =  0 
second  the  equation  of  motion  for  the  current  density 
m  1^  J(Rt)  =  -  V»^(Rt)  -  p(Rt)  VV^(Rt)  -  J  dx  VV(R  -  x)  <n(Rt)  n(xt) 

The  last  term  in  the  above  equation  describes  through  the 
correlation  function  <n(r,t)n(x,t)>^  all  the  complicated  effects  of 
the  Pauli  and  Coulomb  hole  surrounding  each  electron  (in  the  presence 
of  the  external  field).    In  lieu  of  an  exact  evaluation  of  this  term 
or  main  aim  is  to  extract  a  local  field  correction.    Noting  that 

<n(Rt)  n(xt)>^  =  n(Rt)  n(xt)  {g(R,  x;  t)  -  1} 
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where  g(r,x;t)  is  the  non-equilibrium  pair  correlation  function,  the 
simplest  approximation  is  to  replace  this  function  by  its  equilibrium 
value  g(r  -  x),  thus  accounting  for  the  local  depletion  of  charge 
density  but  neglecting  the  dynamics  of  the  hole. 

This  approximation  replaces  the  Hartree  potential  by  a  local 
effective  potential  given  by 

2 

Vg^^(R«)  =  V^(Ko))  -  ^  G(R)  p(Ku) 


wi  th 


G(R)  =  -  i  ;  [S(K  -  q)  -  1]  (5.2) 

Relating  Tr^p(k,w)  to  thas  eaaectave  potentaal  an  the  same  manner 
as  a  free  particle  calculation,  namely  as 

P  P.  f  (p  -  K/2)  -  f  (p  +  K/2) 
pa  to  -  p.^/m  +  ie 

then  the  continuity  and  current  density  equations  together  yield 

xCKco)  =  X-(Kco)/l  -  ^  (1  -G(K))  x^CKco) 
"  0 

where  Xo(k,w)  is  the  Lindhart  polarizabi 1 i tyl53  of  an  ideal  Fermi  gas 

f^(p  -  K/2)  -  f  (p  +  K/2) 

XqCKco)  =  I   :  5  

pa  CO  -  p»^/m  +  ie 

Although  G(k)  contains  the  unknown  structure  factor  S(k)  of  the 
electron  fluid,  this  can  be  determined  by  requiring  consistency  with 
the  structure  factor  obtained  through  the  use  of  the  fluctuation 
dissipation  theorem.    Of  use  numerically  is  the  little  known  direct 
relation  between  G(q)  and  the  pair  distribution  function  g(r)^^'* 
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00 


G(0)=  1  -  0  J   dr  g(r)  j^CQD 


0 


and  by  using  orthogonality  of  the  spherical  bessel  functions 


2r 


2  0° 


g(r)  =  1  - 


;   dO  0  G(0)  j^(Qr) 


0 


These  results  are  extremely  useful  because  they  enable  one  to  generate 
a  G(q)  directly  from  the  pair  correlation  function  and  to  test  G(q)  by 
constructing  the  corresponding  g(r). 

In  the  STLS  scheme  both  exchange  and  correlation  corrections  to 
the  Hartree  field  are  automatically  taken  into  account  through  its 
self-consistency.    Moreover  it  is  straightforward  to  recover  the 
Hubbard^^^  result 


by  substituting  for  the  structure  factor  in  eq.  5.2  its  Hartree-Fock 
value.    The  Random  Phase  Approximation,  which  is  the  dynamic  extension 
of  the  Debye-Huckel  theory,  can  be  viewed  in  this  scheme  as  the 
approximation  that  g(r,x;t)  =  1 ,  and  leads  to  the  result  G(k)  =  0. 

At  zero  temperature  a  further  constraint  can  be  used  to  check  the 
validity  of  any  approximate  microscopic  theory.    The  compressibility 
"sum  rule"  states  that  the  long  wavelength  limit  of  the  static 
response  should  agree  with  the  result  obtained  for  the  compressibility 
from  differentiation  of  the  ground  state  energy.    (This  and  other 
ground  state  properties  of  the  quantal  electron  gas  are  covered  in 
Ichimaru.^    Both  the  Hubbard  and  STLS  equations  violate  this 


G(K)  = 
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constraint.    This  has  led  to  modifications  of  the  Hubbard  correction 
factor'"-'" 

G(K)  =  -ir^  X 

r  +  aKp 

and  of  the  STLS  scheme  by  Vashista  and  Singwi^l 

G(K)  =  (1  +  a  |-)  {-  ^  ;  ^  ^  [S(K  -  Q;  p)  -  1]}  (5.2) 

where  the  additional  free  parameter  "a"  is  to  be  varied  such  that 

compressibility  sum  rule  is  nearly  satisfied.    However  it  has  been 
79 

shown     that  if  the  correction  factor  G(k,w)  is  static,  it  is 
impossible  to  satisfy  the  compressibility  sum  rule  and  the  third 
frequency  sum  rule  simultaneously  (the  f  sum  rule  being  denoted  as  the 
first  moment) . 


The  Relaxation  Function 


As  an  alternative  approach  to  calculating  the  correction  factor 
G(k,w),  which  unlike  the  methods  covered  so  far  easily  generalizes  to 
finite  temperatures,  we  shall  consider  approximations  to  the  Kubo 
function    provided  by  the  Mori  continued  fraction  method. 

In  Chapter  I  the  Kubo  function  was  identified  from  the  difference 
of  the  dynamic  and  static  density  response  functions: 

X(KZ)  =  x(K)  +  ipz  C(KZ) 

In  this  section  we  shall  demonstrate  that  this  follows  from  the 
definition  of  the  Kubo  function  as  the  equilibrium  averaged  density 
auto-correlation  function 
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,    «>         >^H     _  -XH 
C^j^Cr  -  r' ,  t  -  t')  E  ^  ;   dX  <e     n(rt)  e  n(r't')> 

0 

-  <n(rt)>  <n(r't')> 

where  <  >  stands  for  the  thermal  expectation  average. 

It  is  easily  shown  that  the  Kubo  function  vanishes  at  large 
times;  spatial  Fourier  transformability  is  insured  by  the  subtracted 
term  <n><n>.    The  intimate  relationship  of  the  Kubo  function 
with  the  density  correlation  function  can  by  exposed  by  writing  the 
Kubo  function  in  the  form 

C(r  -  ?'.  t  -  t')  =  <n(rt)  n(r't')>  -  <n(rt)>  <n(r't')> 
where  we  defined  the  Kubo  transform  of  an  operator  as 

J  E  1  ;   dX  e^"  J  e"^" 
'  0 

By  taking  matrix  elements  diagonal  in  the  Hamiltonian  of  the  Kubo 
transformed  operator,  it  is  easily  seen  that  the  Kubo  function  reduces 
in  the  high  temperature  limit  to  the  density  correlation  function. 
(The  temperature  =  0  limit  however  is  ill-defined).    Other  notations 
abound  in  the  literature;  in  lieu  of  the  Kubo  transform  we  shall 
employ  the  semicolon  bracket^"^^ 

,    P  XH  -XH 

<A;  B>  E  ^  ;^  dX  Tr  {p^^  e     A  e  B} 

Using  time  translation  invariance  and  the  cyclic  property  of  the 
trace  one  finds  that  the  Kubo  function  has  the  property 

^  ft  "^^^  -  =  p  ^"^"^  - 
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The  time  Fourier  transform  of  this  equation 
c(r  -  r'.  CO)  =  I  X"  ^'^  ~ 

yields  as  a  consequence  that  the  static  response  is  related  to  the 
time  =  0  limit  of  the  Kubo  function 

^(K)  ^  J  d«  ^„  moi  ^  J.      3C(a,)  =  pC(t  =  0) 

As  a  further  consequence  we  see  that,  unlike  the  intimately  related 
density  correlation  function,  the  frequency  moments  of  the  Kubo 
function  are  those  of  the  response  function  x"  without  encumbering 
factors  of  coth(bhw/2). 

Taking  the  Laplace  transform  of  the  time  derivative  equation 
yields  the  relation  referred  to  in  Chapter  I 

X(KZ)  =  x(K)  +  ipz  C(KZ)  (5.3) 

alternatively  this  can  be  viewed  as  simply  the  integration  by  parts  of 
the  Laplace  transform  of  the  Kubo  function  (note  the  tilda) 

00 

C(rZ)  =  ;   dte^^^  C(rt)  In  Z  >  o 

0 

P  d(o_  C((o')       1  p  dco   x"(r(o)         Tn  7  .  n 
-  J  2Tr.  CO'  -  Z  -  3  J  iri  to(co  -  Z)         ^  Z  "  0 

The  last  equality  follows  from  the  integral  representations  of  the 
complex  and  static  response  in  terms  of  the  response  function. 

Based  on  analytic  properties  the  Laplace  transform  of  the  Kubo 
function  can  always  be  represented  as^'^ 
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with  D(z)  an  analytic  function  for  Imaginary  (z)  not  zero.    It  is  the 
function  D(z)  that  we  will  investigate  further  and  which  we  shall  see 
can  be  expressed  in  the  form  of  a  continued  fraction  that  may  be 
truncated  at  some  appropriate  point  of  approximation. 

Scalar  Products  of  Operators 

We  start  by  considering  in  detail  the  time  evolution  of  physical 
observables.    In  quantum  mechanics  operators  satisfy  the  equation  of 
motion 

^  A(t)  =  j:^  [A(t).  H]  E  iLA(t) 

which  defines  a  linear  Liouville  "super-operator,"  that  is  an  operator 

on  operators.    If  we  write  this  out  in  the  form  of  matrix  elements, 

taken  for  a  complete  set  $.  of  quantum  states,  the  first  half  of  the 
above  equation  reads 

f?  AMM(t)  =  |t  (A„  (t)  H     -  H„  a  (t)) 
ot   MN         Tn     mv       vn      mv  vn 

(summation  on  repeated  indices  implied)  or 

If  one  now  treats  the  pairs  mn  or  uv  as  one  index  each,  then  the 
A|jy(t)  can  be  stacked  to  form  a  column  vector,  while  the  Liouville 
super-operator  within  the  braces  of  the  above  equation  can  be  put  into 
the  form  of  a  matrix.    The  point  to  all  of  this  is  that  the  usual 
conception  of  quantum  mechanics  as  a  Hilbert  space  of  say  N  dimensions 
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can  also  be  easily  envisaged  as  a  super-Hi Ibert  space  of  N 
dimensions,  where  the  dynamics  occur  through  a  linear  matrix  operator. 
In  this  super  space  we  can  define  a  bracket  such  that 

(FIG)  E  ^  ;   X  <e  "  F  e    "  G>  =  <F; 
'  0 

Although  often  omitted  the  factor  one  over  beta  is  retained  in  this 
definition  in  order  that  the  bracket  have  the  same  dimensions  as  the 
operator  elements,  and  that  the  bracket  be  well  defined  in  the  beta 
equal  infinity.    The    dagger  representing  Hermitian  conjugate  within 
the  thermal  expectation  brackets  allows  us  to  write  the  spatial 
Fourier  transform  of  the  density-density  Kubo  function  as 

C(Kt)  =  1  (n^(t)  I  n^(o)) 

(noting  that  <n|^(t)><n|^(0)>  is  zero  for  nonzero  wavevector.)  In 
addition  this  scalar  product  satisfies  the  defining  properties  of  a 
Hi Ibert  space: 

(FIG)*  =  (GIF) 

(FIF)  >  0 

(aF  +  bGIH)  =  a(FlH)  +  b(GlH)     a.b  complex  constants 

i.e.  linear  in  bra  (antilinear 
in  ket) 

Furthermore  it  has  the  remarkable  property 


■'1 
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(FIG)  =  (G^'IF"') 

This  formalism  affords  us  the  flexibility  of  manipulations  borrowed 

from  familiar  quantum  mechanics.    For  example  the  auto-correlation 

function  C^^(t)  can  be  written  as 
aa 

Cj^(t)  =  <A(t)  I  A(o)>  =  <e^*'-  A(o)  I  A(o)>  =  <A(o)  I  e"""*""  A(o)> 

=  <A  I  e"^*^  I  A> 

The  Hermiticity  of  the  Liouvillian  follows  from  the  time  translational 

invariance  of  the  brackets.  (Henceforth  we  will  drop  the  argument 

zero,  understanding  A  =  A(t  =  0).) 

Our  aim  is  to  find  an  exact  description  of  the  evolution  of  the 

auto-correlation  function  C,^(t),  in  particular  for  the  case  where  A 

aa 

represents  the  density  operator.    In  achieving  this  a  geometric 
interpretation  of  the  problem  proves  to  be  very  helpful.    In  the  above 
equation  we  can  interpret  the  operator  exp/-itL/  as  rotating  the 
vector  |A>,  and  except  for  a  normalization  constant,  the  correlation 
function  C(t)  is  the  component  of  this  rotated  vector  parallel  to  the 
original  direction.    It  is  therefore  suggestive  to  introduce  a 
projection  operator 

P  =  |A>  <A|A>"^  <A|  =  1  -  Q 

into  the  formal  expression  for  the  Laplace  transform  of  the 
correlation  function 


-84- 


C(Z)  =  <A 
=  <A 


1 


Z  -  L 
i 


A> 


Z  -  LQ  ■*■  Z  -  LQ  "-^  Z  -  L 


L  LP  ^ 


A> 


Here  we  have  used  a  formal  operator  identity  to  reexpress  the  inverse 
of  z  -  LQ.    This  allows  us  to  manipulate  the  Laplace  transform  as 
follows.    The  first  term  of  the  above  equation  yields 


<A 


Z  -  LQ 


A>  =  1  <A  I  A>  =  1  C(t  =  0) 


as  can  be  easily  demonstrated  by  expanding  in  powers  of  1/z  and  using 
the  fact  that  Q|A>  =  0.    As  to  the  second  term  we  insert  the 
definition  of  the  projection  operator 


—  I  A  >  =  I  A  >  T^r-hrr  <  A 


Z  -  L 


<A  I  A> 


~  I  A  >    =  I  A  >  C~^t  =  0)  (Z) 


Z  -  L 


multiply  our  whole  equation  for  the  Laplace  transform  C(z)  by  z  and 
rearrange: 


C(Z)  = 


±3. 


-1 


Z  -  <A 


Z  -  LQ 


A>  C'(t  =  0) 


PC  (t  =  0) 


Noting  that  pC(t  =  0)  =  x  we  finally  arrive  at  the  form  of  eq. 
5.4;  additionally  we  now  have  explicit  knowledge  of  the  unknown 
function  D(z).    This  function  is  reminiscent  of  the  "self  energy"  in 


quantum  mechanics 


<A 


158 


A> 


It  can  be  simplified  further  because 


1 


2  _  LQ  I  -  =  <A  I  {1  +  LO  ^-TIq}  UA> 


=  <A(L)A>  -  i<A 
=  <A  I  A>  -  i<A 


LQ 


OL 


Z  -  OLO 


A> 
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Since  the  second  term  vanishes  as  z  tends  to  infinity,  the  first  term, 
which  we  have  separated  off,  is  the  self-energy  in  the  infinite 
frequency  limit.    In  this  sense  it  bears  a  correspondence  to  the 

1  RQ 

Hartree-Fock  part  of  the  self  energy  in  Green's  function  theory. 
The  second  term  is  referred  to  in  non-equilibrium  statistical 
mechanics  as  the  memory  function  (it  is  nonzero  for  small  z/large 
time) . 

The  fundamental  result  of  this  section  is  that  this  memory 
function  is  itself  of  the  structure  of  a  time  auto-correlation 
function,  namely  of  A.    There  are  two  differences  however:  first, 
only  the  part  QIA>  which  is  orthogonal  to  IA>  enters  in  this 
correlation,  and  second,  the  dynamical  operator  QLQ  is  not  the  full 
Liouville  operator  but  has  projected  from  it  the  part  which  determines 
the  intrinsic  fluctuations  of  the  variable  A.    In  other  words,  if  the 
dynamics  of  the  property  A  is  of  interest  and  we  call  the  many  other 
degrees  of  freedom  the  bath,  then  the  part  OIA>  determines  the 
coupling  of  A  to  the  bath,  and  QLQ  describes  the  internal  dynamics  of 
the  bath  which  feeds  back  to  impose  its  behavior  on  A  itself. 

As  the  memory  function  is  a  form  of  auto-correlation  function 
itself,  we  can  follow  the  same  procedure  as  above  by  projecting  out 
the  time  evolution  of  QIA>  perpendicular  to  its  original  direction. 
This  procedure  can  be  repeated  ad  infinitum  to  obtain  a  continued 
fraction  representation  of  the  memory  function. 

He  should  also  note  that  an  extension  of  the  method  is  to 
generalize  it  to  the  multivariable  case,  where  the  dynamic  property  of 
interest  is  not  a  single  fluctuating  property  of  the  system  but  a  set 
of  independent  observables  A^,  A2,  ...  A^.    When  calculating  the 
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density-density  response  it  is  customary  to  include  the  longitudinal 
component  of  the  current  density  as  well  as  the  number  density.  Such 
an  extension  will  be  considered  later. 

Microscopic  Theory 

We  will  now  use  the  methods  of  the  above  section  to 
microscopically  derive  the  dynamic  response.    Let  us  turn  our 
attention  therefore  to  two  functions,  the  off  diagonal  dynamic 
response  with  momentum  variables  £  and  £' : 

00 

^P'  (dZ)=-l;^dt  e-      <|l^  CnpQ(t).  np.Q]> 
and  the  relaxation  function 

00 

App,(QZ)  E  ;^  dt  e"^*  <npQ(t):  np,Q(o)> 
where  the  off  diagonal  single  particle  density  operator  is 

(For  simplicity  of  notation  we  will  omit  explicit  consideration  of 
spin.)    Of  course  the  actual  dynamic  response  we  are  interested  in  is 
given  by 

and  the  off  diagonal  relaxation  and  response  functions  are  related 
through  a  simple  integration  by  parts,  in  the  same  manner  as  their 
diagonal  counterparts  (see  eq.  5.3) 
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PP" 


Xpp.(OZ)} 


(5.5) 


According  to  the  continued  fraction  approach  presented  in  the 
preceding  section  this  relaxation  function  may  be  expressed  in  the  form 


App,(OZ)  .  (P) 


1 


P")  (P"  I  X(0)  I  P')  (5.6) 


Z  -  io)  +  ?  (Z) 

where  we  have  defined  operators  with  matrix  elements  given  by 
(p  I  i(Q)  I  P')  =<npQ  j  np.Q>EXpp.(Q) 

(p  I  ito  I  p')  E      <hpQ  j  np..Q>  (p"  I  x"^Q)  I  p') 


(5.7) 
(5.8) 


(p  I  ni)  p')  =  1   <ft_n  I  0  — 1—  Q  I  ^.,n>  (P"  I  x-''(0)  I  P') 
p"  '      Z  -  QLQ      '    P  ^ 

We  see  that  the  off  diagonal  momentum  variables  serve  the  utilitarian 
purpose  of  matrix  indices.    As  we  saw  in  the  previous  section  the 
operator  ^{z)  represents  the  back  effects  of  the  bath  on  the  dynamics 
of  n    ;  in  what  amounts  to  neglecting  collisions  between  particles,  we 
will  assume  that  'i'(2)  =  0  (our  first  approximation).    This  allows  us 
to  use  the  identity 

_J          1    r 

-  =  2  I ^  Iw  +  1/ 
Z  -  iu         Z  -  iw 


and  so  combine  eqs.  5.5-5.7  to  yield 


1 


p")  (p"  I  io)  X  I  p') 


(5.9) 


Z  -  lu 

Let  us  consider  first  the  operator  x(q)  or  equivalently  the 
function  Xpp,(q).    This  function  contains  more  information  than 
that  contained  in  the  static  density  response  function  x(q)  supplied 
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for  example,  from  the  integral  equations  of  Chapter  II.    This  is 
easily  seen  by  the  relation 

N  I  I,  Xpp.(O)  =  I  <nQ:  nQ>  =  x(0)  (5.10) 

P  P 

That  is  the  off  diagonal  information  is  integrated  out.    Thus  to  obtain 
an  appropriate  form  for  the  off  diagonal  response  x    , (q)  we  must 
rely  on  its  more  fundamental  definition  in  terms  of  a  temperature 
Green's  function^^^  whose  Dyson's  equation^^^  is  approximated  by 
assuming  the  vertex  function  "F"  is  dependent  on  the  wavevector  q 
only.    As  we  will  show  this  is  correct  in  the  classical  limit  and  is 
approximately  true  in  general  for  fermions.^^^    This  (second) 
approximation  reduces  x    , (q)  to  the  form 

^pp'^^^  =  ^p^^^  ^pp'  -  I^^^Q^  %'^Q^  ^5.11) 

where 

.Ret  ,  ^  ^Ref 


00 


do(0)  =  WZ  !     de  -TT         am  —   ,  , 


and 


pret.  V  -  ]  

p         ~  0  

e  -  Ep  +  |i  -  Ep(e) 

is  the  retarded  single  body  Greens  function,  e°  the  kinetic 

P 

energy  of  a  free  particle  and  Z  (b)  the  self-energy. 

The  function  d^(q)  is  easily  evaluated  for  two  cases.  If  the 
self-energy  z  (e)  is  independent  of  e  we  have 

d  (Q)  =  _  1       P  +  0/2  P  -  Q/2 

P  '        P  +  Q/2  "  ""p  -  0/2 
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where 

-p     -p  •  -p 


c  0  v 


and 

f(e)  =  1/[1  +  el^^"  - 

is  the  Fermi  distribution  (note  that  a  grand  canonical  ensemble  is 
necessary  for  this  portion  of  the  evaluation  only).    In  the  other 
case,  namely  the  classical  limit,  d  (q)  reduces  to  the  normalized 
Maxwel  1-Boltzmann  momentum  distribution  function  (j)|^j^(p): 

lim  1  dp(Q)  =  <t>^g(p) 

and  so  eq.  5.11  agrees  in  the  classical  limit  with 
1  1     N     iOr.  N     iOr  . 

=  •t'MB^P)  «(P  -  P')  +  <I>MB^P^  ^S^Q)  -  1>  *MB^P'^ 

where  <->^  denotes  the  classical  average  over  the  canonical 
ensemble.    This  verifies  the  classical  limit  of  the  vertex  function  as 
being  dependent  on  the  wavevector  q  only  and  identifies  it  as  S(q)-1. 

In  general  the  vertex  function  can  be  identified  by  requiring 
agreement  with  eq.  5.10.    We  obtain 

where  we  have  introduced  a  reference  static  density  response  by 
defining 

XqCQ)  =  I  I  d  (Q) 
P 
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In  the  classical  limit  x^Cq)  =  1  and  x(q)  =  S(q)  so  we  recover 

the  classical  result  eq.  5.12.    In  the  quantal  limit  when  we  neglect 

self-energy  effects,  x^Cq)  is  simply  the  static  limit  of  the  Lindhart 

o 

or  RPA  polarizabi lity. 

Next  let  us  consider  the  matrix  elements  of  the  operator  W  . 
First  we  note  that  the  above  results  for  the  operator  x(q)  allow  us 
to  write 

also  we  note  that  we  can  exactly  evaluate 

<^o'  "p'o"  =  -  U        -  pp.o/2>  ^pp'  =  V^^  W  ^'-^^^ 

in  terms  of 

namely  the  momentum  distribution  of  particles  in  the  interacting 
system.    We  have  introduced  the  function  b  (q)  for  notational 
conveni  ence. 

Inserting  these  results  into  eq.  5.8  yields  the  result 

(P  I  ico  I  p')  =  -  i  flpQ  5pp.  -  n^  CqJQ)  Ibp(O)  (5.15) 
Here  we  have  defined  an  effective  frequency  "omega" 
%Q  -  -  bp(Q)/dp(0) 

which  is  motivated  by  the  fact  that  in  the  limit  of  zero  self-energy, 
where 

Pp  =  f(^p) 
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we  may  write 

'  %  =  i  ^^p+Q/2  -  ^?-Q/2>  =  ?,  P  • 
We  have  also  introduced  a  quantal  direct  correlation  function 

in  agreement  with  the  results  of  Chapter  II.  We  note  again  that  in 
the  appropriate  limit  this  function  reduces  to  the  classical  direct 
correlation  function.    It  now  follows  from  eq.  5.15  that 

(P  I  {z-  i^}-^  I  p')  =  )  - 

pq 

bp(Q)     _J   npCQ^(Q)  ^ 

^    N    >  ^z.iLpQ^  ^1  .  n^QQ,(Q)  x,(Oz)>  (5.17) 

where  the  reference  dynamic  response  defined  as 

1 

is  consistent  with  the  reference  static  response  defined  earlier.  In 
general  x^Cq)  is  a  complicated  functional  of  the  self  energy.  But 
if  we  make  the  assumption  (third  approximation)  that 

p  =  <a  ,  a  >  =  f(e°) 
^p       p'    p      ^  p^ 

then  this  reference  response  reduces  to  the  ideal  fermion  or  Lindhart 
dynamic  response  function.    We  will  see  later  on  that  although  the 
pair  distribution  function  changes  readily  from  its  ideal  fermion  form 
as  interparticle  interactions  are  "turned  on"  the  momentum 
distribution  is  in  fact  an  extremely  insensitive  function,  and  so  this 
assumption  is  in  fact  adequate. 
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By  using  the  results  eqs.  5.17,  5.8,  5.14  and  summing  equation 
over  the  variables         a  few  easy  manipulations  yield  the  result 

X.(Qco) 

We  note  that  the  corresponding  classical  result,  where  the  effective 
potential  has  been  replaced  by  the  direct  correlation  function,  has 

1  CO     1  "30 

been  independently  derived  elsewhere     '      via  a  variational 
criterion  to  find  approximate  eigenfunctions  and  eigenvalues  of  the 
classical  Liouville  operator. 

The  above  equation  is  the  main  result  of  this  chapter.  Through 
it  and  the  fluctuation  dissipation  theorem  we  are  able  to  obtain  the 
pair  distribution  function  from  the  corresponding  analogue  function 
provided  by  the  integral  equations  of  Chapter  III,  and  in  a  manner 
where  all  approximations  are  systematically  accounted  for. 


CHAPTER  VI 
NUMERICALLY  EVALUATING  S(q) 


This  section  details  the  procedures  for  calculating  the  pair 
distribution  function  g(r)  utilizing  the  results  for  the  "analogue" 
pair  distribution  function  n(r|v)/n^  output  from  the  quantal 
integral  equations  discussed  in  Chapter  II. 

We  employ  the  fluctuation  dissipation  theorem 

00 

S(Q)  ^  J.  ^  diMi  _  2  =  CO  +  ie)  (6.1) 

together  with  the  approximation  derived  in  the  last  chapter 

x.(Qz) 

x(Qz)  =   ^   (6.2) 

where  the  quantum  mechanical  direct  correlation  function  can  be 
obtained  from  the  analogue  pair  correlation  function  directly  in 
spatial  Fourier  transform  space  as 

Cq^(O)  =  C^"(0)/Xq(0) 

It  must  be  remembered  that  it  is  the  analogue  pair  correlation  function 
that  one  obtains  from  the  classical  form  of  the  Ornstein-Zernicke 
relation  in  the  QHNC  system  of  equations.    Note  also  that  we  have 
employed  the  notation  of  Chapter  III:    the  density  response  function 
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is  defined  with  an  additional  factor  of  -1 /beta*densi ty  not  found  in 

7Q 

Chapter  II  or  commonplace  in  the  literature.      As  defined  here  the 
ideal  density  response  is  dimensionless 

(Oz)  -  ^  r   H-  0/2)  ~        -  0/2) 

^0^"^^^  -  p  J  (2Tr)2  1^^^  -  ^^^^^  +  Q/2)  -  3e(K  -  Q/2)]  ^^'"^^ 

g  denotes  the  spin  multiplicity,  rho  the  density,  f  the  Fermi  momentum 
distribution  and  E(q)  the  kinetic  energy  of  wavevector  q.  For 
arbitrary  complex  frequencies  "z"  the  ideal  response  has  both  a  real 
and  imaginary  part.    Using  the  identity 

lem     J     =  P  V  ±  i^S(x) 
e-'O  X  +  ie 

we  see  that  along  the  real  frequency  axis  there  is  a  discontinuity; 
approaching  the  real  frequency  axis  from  above  (below),  we  have 

Xq(Q,  z  =  0)  ±  ie)  =  Xq(Q(o)  ±  IXqCQo) 

Under  the  approximation  we  have  derived  for  our  local  field 
correction  factor  (eq.  5.16,  5.18),  we  see  that  the  interacting 
response  shares  this  property  as  well  by  virtue  of  the  fact  that  C 

qm 

is  purely  real.    This  allows  a  simpler  numerical  evaluation  of  the 
fluctuation  dissipation  theorem  as  follows:^^^    The  integral  over 
the  entire  real  frequency  axis  involving  the  imaginary  piece  of  the 
complex  response  can  be  rewritten  as  a  complex  variable  contour 
integration 


c(0)  -  r     d(f>z)   2   Bx(Oz) 

-         2.     ^  _  ^-3.z  2i 
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where  the  contour  c,  is  given  by  Fig.  6.1. 


Imz 


Rez  = 


CO 


Figure  6.1 
Integration  contour  in  complex  plane. 

By  closing  the  loops  in  the  upper  and  lower  half  of  the  complex  plane 
with  an  infinite  radius  arc  (over  which  the  integrand  contributes 
nothing)  we  can  then  proceed  to  deform  the  contour  to  that  of  Fig.  6.2, 


Imz 


Rez 


(0 


Figure  6.2 
Integration  contour  in  complex  plane. 


This  loop  encloses  the  poles  at 


2Trmi 
3f. 


z  =  =^  m  =  0,  +  1 ,  +  2, 
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arising  from  the  exponential  in  the  denominator,  from  which  we 
evaluate  the  contour  integral  by  using  the  Cauchy  Residue  theorem.  We 
obtain 

00 

S(Q)  =     I     x(Q.  fxo  =  ^) 

This  formulation  is  numerically  advantageous  for  two  reasons.  Firstly 
the  relation^^ 

V(Q)  tanh(^)  S(Qto)  =  -  V(Q)  Im  x(Q,z  =  «  +  ie) 

=  ^"^^9  -  e(Q.Z  Jco  +  ie) 

illustrates  the  phenomena  of  plasmon  resonance;  that  is,  for  real 
frequencies  and  small  wavevectors  q  the  dynamic  density  correlation 
function  S(q,w)  behaves  like  a  delta  function  in  frequency  due  to  the 
vanishing  of  the  imaginary  part  of  the  dielectric  function  epsilon. 
The  summation  formula  avoids  those  complications  arising  from 
numerical  quadratures  over  delta  peaked  integrands.    Secondly  along 
the  imaginary  frequency  axis  the  ideal  response  is  purely  real  and  an 
even  function  of  (purely  imaginary)  frequency.    For  finite 
temperatures  this  allows  us  to  write 

00 

S(Q)  =  t{Px(Qo)}  +  2t    I    (5x(Q,  2TrmTi) 

m=l 

where  tau  =  1/beta.    In  the  zero  temperature  (tau  =  0)  limit  the 
product  "3x"  remains  finite,  and  the  sum  is  a  function  sampled  at 
infinitesimally  closely  spaced  points.    Using  Gregory's  formula^^^ 
this  summation  can  be  written  in  terms  of  a  definite  integral  plus 
correction  terms 
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(2x)    I    3x(Q.2TrTi  +  iJ(2TrT))  = 

00 

=  z  J   3x(Q,iv)  dv  +  ^  px(Q.2irxi)  +  0(t^) 
Zk-z  ^ 

There  are  examples^^'^  where  the  tau  =  0  limit  (integrating  through  a 
continuous  sum)  of  this  asymptotic  expansion  does  not  equal  the  result 

1  " 

i  J   &x(Q,  iv)  dv 

0 

due  to  cancellations  involving  the  lower  limit  of  the  quadrature  and 
the  first  correction  term.    However  later  numerical  results  will 
validate  the  assumption  that  (under  the  local  field  approximation  we 
derived)  the  correct  zero  temperature  form  is 

1  °° 

S(0)  =  I  !   3x(0iv)  dv 

The  Lindhard  Function 


Numerical  evaluation  of  the  expressions  that  we  have  just  derived 
for  S(q)  depend  critically  on  the  evaluation  of  the  free  response 
along  the  imaginary  frequency  axis  for  arbitrary  temperatures  and 
densities.    At  this  point  let  us  investigate  further  some  of  the 
mathematical  properties  of  this  function. 

For  convenience  we  introduce  a  frequency  dependent  parameter  nu 
by  the  expression 

.2 

where 
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,  .2  1/2 
^  =  ^mkt  ^ 

is  the  thermal  wavelength  and  w"  is  the  magnitude  of  the  purely 
imaginary  frequency.    Note  that  nu  has  dimensions  of  a  wavevector 
squared.    Upon  carrying  out  the  angular  integrations  in  eq.  6.3,  we 
obtain 

3X  (Ov)  =  l(^)  ^  f  ^  f(K)  ^  in  {Vi^^ioi-^^} 

The  zero  temperature  limit  of  the  right  hand  side  is  an  analytic 
function  obtained  by  integrating  by  parts 

1  /2rn.  3   4ir     Q       .^^  ..  ,  dfv 
^    'K-    P  (2ir)^  ^0 

The  remaining  quadrature  over  the  Lindhard  function 

u,  _  2K  +  4K^0^  -  0^  rVL^IoLLEKQl!, 

0  4  "■to  o  o/ 

^  40^  +  [Q^  -  2K0]'^ 

-  ^  {tan-l  -  tan-'  2^-^}  (6.5, 

is  trivial  at  zero  temperature  where  the  derivative  of  the  Fermi 
distribution  is  simply  a  delta  function  in  wavevector  k  at  the  Fermi 
surface. 

To  evaluate  the  static  (nu  =  0)  free  response  the  integration  by 
parts  keeps  the  integrand  finite;  note  however  there  is  a  logarithmic 
singularity  in  its  slope  at  q  =  2k.    This  makes  the  numerical 
evaluation  at  finite  temperatures  non-trivial,  especially  for  low 
temperatures  where  the  temperature  dependance  is  rather  weak  and  comes 
from  a  narrow  region  around  the  Fermi  surface.    We  further  note  that 
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regarding  summation  over  frequency  values  (which  occurs  in  our 
approximation  for  the  local  field  correction  factor)  the  convergence 
is  very  slow  at  low  temperatures  and  large  wavevectors.    For  this 
reason  several  attempts  to  analytically  evaluate  the  finite 
temperature  free  response  have  been  made. 

Khanna  and  Glyde^^^  and  Gouedard  and  Deutsch^^  have  presented 
an  exact  evaluation  along  the  real  frequency  axis  by  summing  the 
residues  of  the  poles  of  the  Fermi  distribution.    Their  analysis  may 
be  extended  to  arbitrary  complex  frequencies;  however  the  series  is 
slow  to  converge  and  the  zero  temperature  limit  requires  the  summation 
over  a  continuum  of  infinitesimal ly  spaced  poles.    It  is  interesting 
to  note  how  the  logarithmic  singularity  in  the  slope  of  the  Lindhard 
(zero  frequency  free  response)  function^^^ 

2Kp     4Kp  -  Q  +  2Kp  ^ 

"Q~  *  ^"  ^qT-2K^^ 

dissolves  when  finite  temperature  effects  are  included: 

2Kpa     4K^Y  -  0^  +  (0  +  2Kpa)^  ^ 

~n~  +   y   In  -^-^  ^  +  0(t  ) 

AQ""  AK^b^  +  (0  -  2Kpa)'^ 

where 


T  =  KT/cp  Y  =  p/sp  a  =  ^^Y  +  2^Y^  +  ir^T^ 


h        1        1,2  22 

Wasserman  et  al.^^''  presented  an  evaluation  (again  along  the 
real  frequency  axis)  based  on  a  combination  of  the  Sommerfeld 
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expansion  and  a  Laplace  transform  representation  of  the  Fermi 


distribution 


168 


These  results  can  be  generalized  to  arbitrary 


complex  frequencies  by  extending  the  simpler  approach  of  Kanhere  et 
22 

al .      We  write 


3x(0V)  =  l(^)  ? 


1  1 


LGi}^  +  0)  -  G(7r  -  0)] 


2\2'  p  2/  20 


where 


G(0)  =  -  G(-  Q)  =  J    dK  Kf(K)  In 

0 


0  +  2K 
0  -  2K 


Following  Kanhere  et  al .  we  obtain 


G(y)  =  (z^  -  y^/4)  In 


V  +  2z 
y  -  2z 


24.2  0^ 


6z' 


■2z  +  y' 


2  2 
t  F  (^^    -  V  ) 


in  terms  of  the  dimensionl ess  variables  y  =  q/k^,  t  =  l/beta*e^  and 
2 

z   =  beta*u/e|r.    Here  we  have  introduced 


F^(x)  =  ;   dO  In 
0 


0  +  X 
0  -  X 


eQ.  1 


a  real  function  of  a  complex  variable  (note  the  absolute  value 
braces.)    By  expressing  x  in  terms  of  its  real  and  imaginary  pieces, 
integrating  by  parts  and  expanding  the  resulting  logarithm  we  find 
that  this  integral  has  an  exact  representation  in  terms  of  exponential 
integral  functions  of  complex  arguments 


00  N+1 

F  (x)  =  Real  {    I    ^        e"^^  E,(NZ)  -  e"^^*  F,(-  NZ*)} 

N=l  '  ' 


,169 


The  approach  of  Dharma-Wardana      is  to  note  that  the  Fermi 
function  can  be  written  as  an  integral  over  a  dummy  chemical  potential 
variable 
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f(e,  11.  T)  =  ;    0(v  -  e)  A(ni,  v,  T)  dv 

0 

where  theta  is  the  step  function  and 

A(vi,  V.  T)  =  {2T[1  +  cosh(>i  -  v)/T]}"^ 

is  the  temperature  weight  function.    Inserting  this  representation  of 
the  Fermi  function  into  the  integral  definition  of  the  free  response 
and  exchanging  the  order  of  integration  results  in  a  quadrature  over 
the  dummy  chemical  potential  of  the  temperature  weight  function  and 
the  analytic  zero  temperature  free  response  evaluated  at  the  dummy 
chemical  potential. 

Although  a  systematic  comparison  of  the  various  methods  was  not 
undertaken,  with  the  machine  word  size  and  speed  of  the  Cray  machine 
available  it  was  found  that  a  numerical  evaluation  of  the  integration 
by  parts  formula  (eqs.  6.4  -  6.5)  was  adequate  when  supplemented  by 
certain  asymptotic  expansions.    For  the  nonzero  frequency  case  several 
asymptotic  expansions  can  be  constructed  from  the  Lindhard  kernal  (eq. 
6.5)  depending  on  the  relative  magnitudes  of  the  variables  k,  q,  nu. 
The  expansion  we  found  most  accurate  was 

o^,/^w^      l/'2mv  /I      1     «K^»  «K^»  ,  ,^ 

Px(Ov)  ~  {-  -  ^      „     +  ^  — r—  +  . . .  }  (6.6) 

^  K-      "     of-  «K°»      3a-^  «K°» 

for  large 

Here  the  double  brackets  «k"»  denote  the  spherical  averaged 
moments  of  the  Fermi  distribution  function.    For  the  zero  frequency 
case  we  found  that  for  large  wavevectors 


-102- 


2 

(5X(Q,  0)  ~  +  — 7   —  +  ...  } 

whi le  for  smal 1  q  values 

Px(Q,  0)  =  2^^^  { — - —  +  — - —  +  . . .  } 
^  fi      «K°»       ''^  «K°» 

We  note  that  the  finite  temperature  moments  «k"»  can  be  given  either 
by  a  Sommerfeld  expansion^ ''^  or  by  exact  analytic  expressions  in  terms 
of  hypergeometric  functions. 

Graphs  of  the  free  response  function  for  differing  wavevector, 
frequency  and  temperature  at  a  typical  metallic  density  characterized 
by  r^  =  2  are  presented  in  the  following  pages.    The  first  figure 
set  presents  a  plot  versus  wavevector  q,  where  the  solid  line  is  for 
the  zero  frequency  response  and  the  dashed  curves  are  for  increasing 
frequency.    We  note  that  zero  frequency  is  a  special  case  with 
differing  limiting  form  near  the  wavevector  origin.    Next  we  present 
the  free  response  as  a  function  of  frequency  for  a  few  selected 
wavevectors.    The  dashed  line  represents  the  leading  term  in  the 
asymptotic  expansion  of  eq.  6.6.    We  note  that  for  larger  wavevectors 
the  free  response  dies  away  in  frequency  very  slowly  and  the 
asymptotic  expansion  is  increasingly  more  representative.    Finally  to 
show  the  effects  of  finite  temperatures  we  plot  the  zero  frequency 
response  versus  wavevector  at  several  temperatures  keeping  density 
constant.    The  feature  to  keep  in  mind  is  the  disappearance  of  the 
logarithmic  singularity  in  the  slope.    These  plots  are  followed  by 
graphs  of  the  free  response  as  a  function  of  frequency  for  a  typical 
wavevector  at  various  temperatures.    We  note  that  the  asymptotic 
expansion  is  most  representative  at  the  lower  temperatures. 
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Convergence  Acceleration 

Because  the  free  response  decays  slowly  with  frequency  at  the 
larger  wavevectors,  two  procedures  were  found  necessary  to  produce 
accurate  results  in  computing  S(q).    The  first  procedure  is  to 
subtract  off  the  corresponding  formulas  for  the  noninteracting 
distribution  S^^(q).    This  results  in 

S(0)  -  S^"(0)  =  -  ;  — —  dv 

at  zero  temperature  while  at  finite  temperatures 

(KT)  pC^"(Q)  3x^(Q0)  -   fixAQO)  ^^o^^^m^ 

S(Q)  -  S^"(0)  =  ^         +  (2  KT)    I  -^-r^^  

1  -  PC^"(Q)  n,=  l  pC^"^Q> 

PXq(QO)  '^^o^^^m^ 

where  v^^  represents  the  temperature  spaced  frequency  points  along 
the  imaginary  axis.    The  noninteracting  distribution  can  be  accurately 
computed  separately  and  the  difference  now  decays  in  frequency  as  the 
free  response  squared.    We  next  note  that  for  large  q  and  or 
frequencies  the  denominators  in  these  expressions  quickly  reach  a 
value  of  unity.    For  this  reason  we  can  accelerate  convergence  by 

Y~  -p 

subtracting  a  reference  S     (q)  defined  at  zero  temperature  as 

cref,^.     pC^"(0)  1      r.^REF.^    .,2  . 
^  =  3Xq(Q0)  ^^'^^^ 

while  at  finite  temperature 
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The  reference  response  function  in  these  formulas  is  any  expression 
that  closely  approximates  the  actual  free  response  and  allows  an  exact 
evaluation  of  the  summation/integration  in  frequency  appearing  in  the 
above  formulas. 

At  zero  temperature  the  reference  response  can  be  taken  as  the 
exact  closed  form  expression  for  the  free  response  itself.  The 
integral  over  frequency  of  its  square  is  a  universal  function  of 
q/k.^  (k^  the  Fermi  wavevector).    However  it  is  not  known  in  a 
closed  form  expression.    This  universal  function  was  numerically 
evaluated  and  fitted  piecewise  by  polynomial  approximations.  The 
large  q  tail  was  found  to  be  accurately  represented  by  the  first  two 
terms  in  an  asymptotic  expansion  of  the  form  a/q"2+b/q"4+. . . 

At  finite  temperatures  we  do  not  have  a  closed  form  expression  to 
work  with  so  we  used  the  first  term  in  the  asymptotic  expansion  of 
eq.  6: 

A     Q   +  ("2)    ^l^^^ni         z   +  m 
A 

where 

This  allows  the  summation  over  the  reference  response  function  squared 
to  be  written  at  any  temperature  as 

I  [px'^^^O  n  )]2  =  ^  ^(Z) 
m=l  f,^ 
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where 

<J>(x)  =  I      2^2  2 
m=l  (x'^  +  m^r 

is  a  universal  function.    This  function  has  an  exact  summation.  It 
can  be  obtained  from  the  partial  fraction  expansions  of  cothz  and 
cothz  (use  the  fact  that  cschz  =  icsc(iz)).    We  have 

(cschz)    =  "2  +  2  {  — 2  2  — 2  2~2  "*■' 

Z  (Z+ir)  (Z+4ir) 

1  1  72       2  72      y,  2 

2  =  3  +  2{-2  +  ._2      ,  2,2  ^  ••'> 

Z  (Z+ir)  (Z+4ir) 

00 

and    (cschz)^  +  ^  cothz  =  21^   I  — = — ^  5-r 

^  m=-oo  [Z^  +  (mrr 

An  analytic  expansion  about  the  origin  can  be  written  in  terms  of  the 
Riemann  zeta  function  of  integer  argument  zeta^ 

<t>(x)  =   +  (5   -  1)  -  2x^(5.  -  1)  +  3x^(^„  -  1)  -.-  ... 

Its  asymptotic  form  can  be  obtained  by  using  Gregory's  formula^ to 
express  the  sum  as  an  integral  plus  correction  terms.    By  itself  it 
leads  to  a  bad  representation,  but  by  expanding  all  the  correction 
terms  in  powers  of  x  and  regrouping  one  obtains  a  very  accurate 
asymptotic  expansion 

Ax^     2x^  x^ 

We  note  in  passing  that  using  the  approximate  form  of  eq.  7  for  the 
reference  response  at  zero  temperature  results  in 
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00 

which  clearly  indicates  the  utility  of  subtracting  off  the  reference 
S(q)  breaks  down  for  small  wavevectors. 
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Figure  6.3 

One  over  temperature  times  the  dimensionless  free  response 
function  [1/ryd]  is  plotted  as  a  function  of  wavevector  q  [1/bohr]  in 
the  limit  of  zero  temperature  and  density  of      =  2.    The  solid  line 
represents  zero  frequency;  the  dashed  lines  represent  (imaginary) 
frequencies  of  1.0,  0.5  and  0.25,  0.10,    in  rydbergs. 
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Figure  6.4 

One  over  temperature  times  the  dimensionless  free  response 
function  [1/ryd]  is  plotted  as  a  function  of  (imaginary)  frequency 
[ryd]  in  the  limit  of  zero  temperature  and  density  of      =  2.  The 
solid  line  represents  values  for  wavevectors  of  1.0,  1.5,  2.0,  and  2.5 
[1/bohr]  while  the  dashed  line  represents  an  approximation  based  on  an 
asymptotic  expansion. 
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Figure  6.5 

One  over  temperature  times  the  dimension! ess  free  response 
function  [1/ryd]  is  plotted  as  a  function  of  wavevector  [1/bohr]  at 
zero  frequency  and  density  rc  =  2.    The  plots  show  the  variation 
with  temperatures  of  0.0,  0.25,  0.50  and  1.0  in  units  of  Fermi  energy. 


Figure  6.6 

One  over  temperature  times  the  dimensionl ess  free  response 
function  [1/ryd]  is  plotted  as  a  function  of  (imaginary)  frequency 
[ryd]  at  a  wavevector  q  =  1.0  [1/bohr]  and  density      =  2.  The 
plots  show  the  variation  with  temperatures  of  0.0,  0.5,  1.00  and  2.00 
in  units  of  Fermi  energy.    The  dashed  line  represents  an  asymptotic 
approximation  which  is  independent  of  temperature. 


CHAPTER  VII 
THERMODYNAMIC  PROPERTIES 


In  this  chapter  we  will  present  results  for  g(r)  based  on  the 
QHNC  QOZ  equations  and  the  approximation  for  the  local  field 
correction  factor  derived  in  the  last  chapter.    These  will  be  compared 
with  the  results  obtained  under  the  binary  Slater  Sum  approximation 
and  those  of  Pokrant. 

In  classical  statistical  mechanics  knowledge  of  the  pair 
distribution  function  g(r)  alone  is  sufficient  for  the  calculation  of 
most  desired  thermodynamic  quantities, ^^'^^^  for  example  the 
internal  energy  per  particle 

e  =  I  kT  +  I  ;  dr  g(r)  V(r) 
the  pressure 

^  =  1  -        ;  df  g(r)  V(r) 
or  the  compressibility 

=  1  +  p  ;  dr  [g(r)  -  1]  =  1  -  p  J  dr  C(r) 


'  3p 


For  the  one  component  plasma  g(r)  must  be  replaced  by  g(r)  -  1  in 
all  formulas  whenever  g(r)  multiplies  the  interaction  potential  v(r). 
This  is  because  the  OCP  is  actually  a  two  component  system  in  which 
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the  second  component  is  taken  to  be  completely  structureless.  It  has 
a  state-dependent  potential  energy 

uCr^  =1         -  1  N  p  J  dr  !^ 

The  actual  potential  v(r)  entering  the  thermodynamic  formulas  is  that 
of  a  Coulomb  potential  minus  its  zero  wavevector  spatial  Fourier 
component  due  to  the  neutralizing  background,  which  is  equivalent 
under  spatial  integrals  to  h(r)  =  g(r)  -  1  times  a  pure  coulombic 
form.    Furthermore  the  compressibility  must  be  calculated  in  terms  of 
the  pair  correlation  function  in  which  c(r)  must  be  replaced  by 
c(r)  +  bv(r)  to  insure  finiteness.    This  can  be  derived  from  a 
microscopic  field  fluctuation  calculation.^^    The  normal  derivation 
of  the  compressibility  in  terms  of  the  fluctuation  in  number  density 
of  a  grand  canonical  ensemble    is  invalidated  by  virtue  of  the  fact 
that  the  potential  energy  in  the  grand  partition  function  is 
state-dependent. ^^^"^^^ 

The  quantal  extension  of  the  mechanical  equations  can  be  obtained 
by  scaling  methods  similar  to  those  used  to  derive  the  classical 
pressure  and  internal  energy  equations^^"^^  but  that  now  require  the 
off-diagonal  density  matrix.    For  example 

PVol  =  ^  lim_^  ;  dr  p^  p(f  |  7. )  _  1  p2  j  ^-g^^^  ^(r) 

where  p(hat)  is  the  quantum  momentum  operator.    On  the  other  hand  the 
compressibility  equation  has  a  purely  statistical  derivation  that 
should  remain  equally  valid  in  quantum  mechanics.    As  the  pair 
correlation  function  asymptotically  reduces  to  the  effective 
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potential       the  criterion  of  finiteness  illustrates  that  the 
analogue  pair  correlation  of  the  QHNC  equation  should  not  be  used  to 
calculate  the  compressibility.    Rather  the  proper  quantal  extension  of 
the  pair  correlation  function  would  be 

C(r)  =  F-^{1/Xq(0)  F  {C^"(r)}}  ~  ^ 

in  agreement  with  the  development  presented  in  Chapter  II. 

The  deviation  from  the  classical  formulas  becomes  apparent  if  we 
write  the  mechanical  state  equations  in  the  compact  form 

E  =  Ejp,,L  '^^'^  <V>^ 
and  from  the  quantum  mechanical  virial  theorem64 
|(PVOL  -  PVol^^^^"-)  =  AK  -  X  <  H  >^ 

Here  the  brackets  subscripted  by  lambda  denote  the  thermal  average 
over  a  hamiltonian  at  intermediate  coupling  H  =  T  +  XV;  H  is  the 
virial  operator  (equal  to  minus  one  half  V  for  coulomb  interactions) 
and 

AK  =  <T>.  -  <T> 
X  0 

is  the  excess  kinetic  energy.    The  excess  kinetic  energy  is  a  measure 
of  the  deviation  of  the  momentum  distribution  of  particles  from  the 
pure  Fermi-Dirac  distribution.    In  principle  it  can  be  obtained  from 
the  pair  distribution  function  computed  at  different  coupling 
strengths  as 

f  =  f  J%X  (3  1^  -  ^)  J  d?  g(f ,  3.  X)  vCr) 
this  expression  follows  from  combining  the  Maxwell  relation 
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with  the  quantum  statistical  definition  of  the  free  energy 
g-|5F(X)  ^  ^    g-3(T  +  \V) 

and  is  equivalent  to 

X 

F(X)  =  F(X  =  0)  +  ;   dX  <v>^ 

0  ^ 

In  classical  statistical  mechanics  the  functional  dependence  of  the 
pair  distribution  function  upon  the  temperature  and  coupling  strength 
enters  through  the  single  (OCP)  parameter 

r  =  pe^/a 

and  so  the  excess  kinetic  energy  vanishes  identically.    In  this  way 
the  excess  kinetic  energy  also  serves  as  a  measure  of  the  quantum 
effects  on  the  pair  distribution  function. 

It  is  important  to  note  that  in  the  QHNC  theory  the  momentum 
distribution  is  not  directly  obtainable— in  principle  we  must 
integrate  over  the  coupling  constant  in  order  to  obtain  thermodynamic 
properties. 

It  is  fortuitous  that  in  the  classical  HNC-OZ  system  of  equations 
the  integration  over  the  coupling  ratio  can  be  done  analytically  and 
the  result  is  expressed  solely  in  terms  of  fully-coupled  distribution 
functions. This  cannot  be  carried  over  into  the  quantal  case 
because  the  effective  potential  depends  on  temperature  and  coupling  in 
a  nontrivial  manner  and  the  QHNC-QOZ  equations  are  for  the  analogue 
pair  distribution;  the  coupling  constant  integration  must  be  performed 
over  the  physical  pair  distribution. 
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At  zero  temperature  the  algorithm  for  the  excess  kinetic  energy 
is  considerably  simpler.    It  is  equivalent  to  an  integration  over 
different  density  values  at  full  coupling: 

£  ^  2^  _  ^       1    J.  s  ^     V(rJ  +  .916}  dr^  [Ryd] 
r  s      r      0  ^ 

This  formula  is  easily  enough  derived  from  the  relations 

T  -  -  g?-      E(r^)       y-'rJr  rlE  (r^) 
s  s  s 

which  can  be  obtained  by  combining  the  definition  E  =  T  +  V  and  the 
virial  theorem 

2  T  +  V  =  -  r  ^ 


s  dr, 
s 

(the  right  hand  side  is  just  the  virial  3*pressure*volume)  while  noting 

"  ■  -  dvSE     -  ^> 

In  addition  at  zero  temperature  we  have  two  more  constraints  available 
to  us.  Since  the  minimum  value  of  T  is  the  noninteracting  Fermi  limit 
we  have 

_3_        rx       2.21  ^  ^ 
-  9     (r^E)  -         >  0 

Ferrel^^^  has  a  more  restrictive  condition 

^  r  V  <  0 
ar^    s  - 

In  principle  the  procedure  is  to  compute  accurately  values  of  the 
interaction  energy  per  particle 
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Pot  E  ^  p  J  dr  ^  [g(r)  -  1] 


at  varying  densities  and  temperatures,  following  either  the  zero  or 

finite  temperature  algorithms  for  calculating  the  excess  kinetic 

energy  as  a  means  of  obtaining  thermodynamic  quantities.  However, 

there  is  an  often  overlooked  drawback.    The  coupling  constant 

algorithm  for  computing  the  energy  of  an  interacting  system  can  be 
1 78 

shown      to  be  rather  unreliable  when  it  is  used  with  approximate 
wavefunctions  or  pair  correlation  functions.    Thus,  among  other 
practical  considerations  such  as  machine  time,  an  alternative  theory 
for  approximately  calculating  the  interacting  momentum  distribution 
might  be  preferable.    In  the  next  section  we  will  compare  results  for 
the  interaction  energy  per  particle  obtained  from  Pokrant,^^  the 
binary  Slater  Sum,  and  QHNC  methods.    Values  for  the  excess  kinetic 
energy  obtained  under  the  BSS  or  Pokrant  methods  will  be  compared  to 
establish  the  magnitude  of  their  contribution  to  thermodynamic 
quantities  and  their  amenability  to  approximation. 


In  either  the  BSS  approximation  of  the  the  theory  of  Pokrant  the 
thermodynamics  are  derived  from  the  classical  formulas  involving  a 
state-dependent  configuration  potential  energy 


Thermodynamics  with  Slater  Type  Effective  Potentials 


e 


-  |3U^(r'\  (i,  p) 
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This  is  just  the  diagonal  Slater  sum  computed  from  a  set  of  properly 
symmetrized  plane  wave  basis  functions  (labeled  by  the  3N  dimensional 
wavevectors  k")  where  the  "A"  is  an  operator  anti symmetrizing  the 
position  coordinates  of  the  spatial  function  to  its  right.  (We 
employed  the  commutivity  of  the  anti symmetrizer  with  the  hamiltonian 
and  its  idempotency  property.)    The  expression  for  the  free  energy 
follows  as 

.Tre-»".  (x3Nk!)-'  ;d-re"^"T 

while  the  Maxwell  relation  for  the  internal  energy  yields,  under  the 
assumption  that  the  effective  potential  is  pairwise  additive,  is 

^  =      +  f  ;  dr  g(r)  ^  [3U^(r.  (J.  p)] 

To  properly  compute  this  integral  we  should  now  split  the  pairwise 
additive  interaction  potential  terms  as 

The  first  piece,  or  direct  part,  arises  from  the  non-commuti vi ty  of 
the  kinetic  and  potential  energy  operators  in  the  exponentiated 
Hamiltonian  sandwiched  between  simple  product-ordered  plane  wave 
states  (arising  from  the  unity  operator  in  a  permutation  expansion  of 
the  anti symmetrizer  operator).    Because  the  potential  energy  operator 
for  the  OCP  is  just  a  coulombic  minus  its  zero  wavevector  component  of 
the  spatial  Fourier  transform,  the  direct  part  of  the  effective 
potential  is  minus  this  component  also.    (This  property  is  invariant 
under  the  multiple  commutations  with  the  momentum  operator.)    As  such 
the  direct  part  has  a  1/r  tail  (the  contribution  of  the  classical 
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limit)  and  the  integral  over  the  product  of  g(r)  and  the  direct 

potential  (minus  a  Fourier  component)  should  be  replaced  by  g(r)  -  1 

times  the  real  space  form  of  the  direct  potential  (i.e.  including  all 

Fourier  components). 

The  second  piece,  or  exchange  part,  arises  from  the  exponentiated 

Hamiltonian  sandwiched  between  simple  product  plane  wave  states 

ordered  in  index  k"  but  with  the  spatial  coordinated  r"  permuted 

out  of  order.    This  results  in  a  function  of  the  r"  that  must  have  a 

zero  wavevector  spatial  Fourier  component,  regardless  of  the  form  of 

the  potential,  as  is  obviously  the  case  when  the  potential  is 
1 79 

absent.        As  such  g(r)  should  not  be  replaced  by  g(r)  -  1. 
The  BSS  results  were  computed  from  the  formula 

M-fjd^  {[g(r)-l]|pU^.g(r)|3U^} 

where  the  direct  effective  potential  came  from  a  calculation  of  two 
antiparallel  spin  electrons  while  the  exchange  effective  potential 
arose  from  parallel  spin  electrons.    A  similar  correction  should  be 
made  to  the  results  of  Pokrant;  however,  the  resolution  of  his 
effective  potential  into  direct  and  exchange  portions  is  tacitly  more 
complicated.    We  note  in  passing  that  the  same  correction  appears  in 
the  pressure  formula 

^  =  1  -  t  ;  df  {[g(r)  -  n  r  1^  3Up  .  g(r)  pu^} 

Interaction  Energy  Results 

Table  7.1  lists  the  interaction  energy  per  particle  obtained  by 
Pokrant,  the  BSS  effective  potential  using  the  above  refinement,  and 
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from  the  QHNC  local  field  correction  factor.    For  comparison  the  same 
quantity  is  also  calculated  using  the  analogue  pair  distribution 
function  obtained  directly  from  the  QHNC  equations;  that  is,  upon 
neglecting  recoil  and  exchange  effects.    That  the  effective  potential 
of  Pokrant  reduces  to  the  BSS  effective  potential  is  reflected  in  the 
high  temperature  behaviors.    At  zero  temperature  the  results  of 
Pokrant  are  in  fair  agreement  with  those  of  Ceperley  and  Alder, as 

AO 

IS  easily  verified  using  the  parametrization. 

r     ZJA      ^916   .2846   ,  ^-i 

E  =     p    -  ^r-  -   :   [ryd]    (for  r   >  1) 

r^^        's       1+1 .0529  vr^  +  .3334  r^  ^ 

A  short  comparison  of  quantum  Monte  Carlo  results  with  those  of 

Pokrant  is  presented  in  Table  7.2.    For  several  r    values  the 

s 

agreement  for  the  interaction  energy  per  particle  is  within  one  half 
percent.    Thus  the  deviation  of  the  BSS  data  from  that  of  Pokrant 
illustrates  the  breakdown  of  the  binary  approximation  at  low 
temperatures.    The  data  calculated  from  the  analogue  pair  distribution 
function  lie  below  the  results  of  Pokrant  at  high  temperatures  but  at 
low  temperatures  approach  those  of  the  BSS  approximation.  This 
clearly  illustrates  the  error  in  approximating  g(r)  by  g^"(r)  as  is 
pervasively  done  in  plasma  theory. However  the  physical  pair 
distribution  function  obtained  from  the  local  field  correction  factor 
calculation  does  not  yield  an  improvement;  the  results  are  uniformly 
higher  than  Pokrant  and  there  is  an  unusual  dip  around  the  temperature 
origin. 

These  results,  although  disappointing,  should  hardly  come  as  a 
surprise  because,  up  to  this  point,  we  have  concerned  ourselves  solely 
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with  the  most  rudimentary  of  the  quantal  extensions  of  the  classical 
fluid  equations.    The  QHNC  equation  does  not  incorporate  the  most 
important  physical  effect  of  many  particle  systems  interacting  with 
long  range  forces,  to  wit  the  screening  of  the  effective  potential. 
The  results  of  more  refined  quantal  integral  equations  will  be 
presented  in  the  next  two  chapters  in  lieu  of  laboring  the  results  of 
the  QHNC  equation  here.    Instead  we  will  consider  more  closely  the 
thermal  properties  of  the  excess  kinetic  energy. 

Kinetic  Energy  Results 

A  comparison  of  BSS  and  Pokrant  values  for  the  excess  kinetic 
energy  is  presented  in  Table  7.3.    Although  the  kinetic  energy 
provides  a  more  stringent  test  than  the  correlation  energy  (by  virtue 
of  the  variational  principle  for  the  ground  state  energy),  it  is 
nevertheless  surprising  that  the  zero  temperature  results  of  Pokrant 
differ  from  those  of  the  quantum  Monte  Carlo  by  fifteen  per  cent  or 
more.    (See  Table  7.2.)    The  total  energies  indeed  are  in  fair 
agreement  (one  or  two  percent)  but  this  is  because  the  excess  k.e.  is 
a  fraction  of  the  magnitude  of  the  interaction  energies  and  because 
the  errors  in  the  two  quantities  are  self-correcting.    This  shows  also 
that  using  the  virial  theorem  to  produce  pressure  values  yields 
inaccurate  results.    (One  must  use  a  method  which  takes  advantage  of 
the  relative  accuracy  of  the  total  energy— i.e.  integrating  energy 
over  inverse  temperature  to  obtain  the  Helmholtz  free  energy  and  then 
differentiating  with  respect  to  volume.)    For  much  the  same  reasons 
the  internal  energies  of  the  BSS  approximation  reproduce  the  results 


-121- 


of  Pokrant  down  to  surprisingly  low  temperatures;  however, 
individually  the  interaction  and  excess  kinetic  energies  differ 
considerably. 

In  spite  of  the  vagaries,  features  of  the  temperature-dependence 
of  Pokrants  excess  kinetic  energy  (which  are  the  only  finite 
temperature  results  available)  should  be  accurate. 

The  feature  that  the  excess  kinetic  energy  is  positive  at  zero 
temperature  is  elementarily  explained  because 

2 

AK  «  J    dp  p^[f(p)  -  f  (p)] 
0 

and  fgCp)  is  a  step  function.    Statistics  constrains  that  any  change 

in  the  distribution  must  be  accomplished  by  an  increased  occupation  at 

a  higher  momentum  value  at  the  expense  of  a  decreased  population  at  a 

lower  value.    A  typical  plot  of  the  interacting  momentum  distribution 

is  presented  in  Fig.  7.1.  °     The  existence  of  a  discontinuity  at 

the  Fermi  surface  (Migdal's  theorem^^^"^^^)  is  modeled  by  the  small 

r^  RPA  theory  of  Daniel  and  Vosko,^^^  extended  by  Gel  dart 

et  al.,^^^  and  at  large  r^  by  Belyakov.^^°   We  note  that  this 

discontinuity  disappears  only  when  the  momentum  distribution  is 

constructed  out  of  orthogonal ized  Wigner  orbitals  at  the  high  r 
191-192  ^ 

Mmi  t. 

At  finite  temperatures  the  momentum  distribution  can  in  principle 
be  calculated  from  the  one-particle  Green's  function^^^  given  the 
self-energy  (related  to  the  two-particle  Green's  function).    At  zero 
temperature  the  self-energy  is  closely  approximated  by  the 
Hartree-Fock  self  energy  plus  a  self-energy  due  to  single  particle 
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coupling  to  density  fluctuations,  the  main  contribution  of  which 

arises  from  coupling  to  the  plasmon  field  J  ^'^"^^^ 

To  model  temperature  effects  we  start  with  the  finite  temperature 

evaluation  of  the  Hartree-Fock  self-energy  following  the  steps  of 

1 59 

Kadanoff  and  Baym.        We  include  spin  degrees  of  freedom  in  the 
equation  of  motion  for  the  single-particle  Green's  function 

2  2 


6(1  T)  -  i  I  ;  drg  V^^(r^  -        g}^\i  2  1'  2+) 
o 


Here  the  second  term  on  the  right  hand  side  is  equivalent  to  the 
self-energy  as  defined  in  terms  of  the  two-particle  Green's  function, 
the  arrows  denote  spin  labels  (a  parallel  equation  exists  for  the 
spin-down  Green's  function)  and  the  shorthand  argument  "1"  stands  for 
the  space-time  argument  of  particle  one,  etc.    The  Hartree-Fock 
approximation  consists  of  decoupling  the  equations  of  motion  by 
assuming 

g{^^  =  G^d  1')  G^(2  2'') 
as  particles  with  opposite  spin  are  distinguishable  and 
g{^^  =  G^d  1')  G^(2  2"')  -  G^d  2"")  G^(2  1') 

with  a  similar  equation  for  the  down-spins.    Because  the  form  of  the 
potential  is  independent  of  spins  we  write 
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^2  ^2 


6(1-T) 


t2  =  fi 


<r^lVg^^|r2>  =  6(r^  -        J        V(r^  -  r3)[-i  G^(3  3)  -  iG^(3  3)] 


-  V(r^  -  r^)  [-  iG^Cl  2)] 


^2  =  ^1 


Now  the  density  distribution  is  given  by  the  diagonal  limit  of  the 
retarded  Green's  function 

n^(r  t)  =  -  i  G^(rt,  rt) 

and  is  a  constant  for  uniform  systems.    As  a  result  the  spatial 
Fourier  transform  of  the  Green's  function  satisfies  the  ideal  Green's 
function  differential  equation 

[ifi  It  -  ^T^P^^  ^T^P'       ^'^  =  - 
where  the  energy  eigenvalues 
_2  , 

^T^P^  =  L  ^  ^T^P^  (7.1) 
are  given  in  terms  of  the  HF  self-energy 

-.1 

E^(p)  =  (n^  +  n^)  J  dr3  ^ir^)  -  ;  V(p  -  p')  n^(p)  (7.2) 

V(p)  is  the  momentum  (spatial  Fourier)  transform  of  the  interaction 
potential.    It  is  the  only  surviving  term  due  to  the  presence  of  a 
neutralizing  charred  background  for  the  OCP.    The  differential 
equation  for  the  ideal  Green's  function  implies  that  the  momentum 
distribution  is  of  the  form 
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3(E.(p)  -  p) 

n^(p)  =  1/{1  +  e     ^  }  (7.3) 

where  the  chemical  potential  is  fixed  by  the  total  particle  number 
densi  ty. 

The  above  eqs.  7.1-7.3  are  to  be  solved  sel f-consi stently  for  the 

finite  temperature  momentum  distribution.    However,  near  zero 

temperature  the  convergence  can  be  quite  slow.    This  is  evidenced  by 

the  fact  that  at  zero  temperature  f(p)  has  a  step  discontinuity  which 

is  hard  to  achieve  through  a  Picard  iteration,  while  the  eigenvalues 

reduce  to  the  Hartree-Fock  ground  state  energy  which  diverge  at  small 
28 

r^.      A  near  zero  temperature  approximation  can  be  derived 
analytically^^^  but  it  is  interesting  to  note  that  the  same  finite 
temperature  system  of  equations  can  be  obtained  from  a  free  energy 
variational  principle:    the  value  F  =  E  -  TS  is  a  minimum  at  constant 
temperature  under  the  constraint  that  the  total  number  of  particles  is 
fixed.    Aside  from  the  kinetic  energy  contribution  to  the  internal 
energy,  the  entropy  can  also  be  expressed  in  terms  of  the  quantum 
mechanical  one-body  Wigner  distribution''^^  as 

— ♦ 

TS  =  -  kT  J  {f  In  f  +  (1  -  f)  In  (1  -  f)} 

(2Tr)'' 

(which  for  spatially  uniform  systems  is  the  momentum  distribution). 
The  potential  energy  depends  on  the  two-body  Wigner  distribution 
functions,  which  however  can  be  split  into  anti symmetrized  products  of 
one-body  Wigner  functions,  incorporating  exchange  effects,  and 
irreducible  two-body  operators.    For  the  OCP  the  exchange  piece  can  be 
calculated  in  terms  of  the  momentum  distribution  in  the  closed  form 
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(2Tr) 


dK  4Tre 


2 


K  +  0 
K  -  Q 


2 

)    f(k)  f(0) 


In  ( 


Retaining  only  this  contribution  and  varying  the  free  energy  with 

respect  to  the  momentum  distribution  yields  the  system  of  equations 

derived  by  the  Green's  function  method.    However  in  this  approach  we 

may  take  an  approximate  form  for  the  momentum  distribution,  dependent 

on  a  number  of  adjustable  parameters  which  we  vary  to  minimize 
1 96-1 97 

F.  This  allows  for  a  fast  accurate  solution  near  absolute 

zero  for  the  excess  kinetic  energy. 

The  result  of  the  temperature-dependent  HF  approximation  is 
presented  in  Table  7.4.    The  data  differ  markedly  from  that  obtained 
by  the  method  of  Pokrant.    The  excess  kinetic  energy  properly  tends  to 
vanish  in  the  classical  limit  under  both  approximations,  but  the 
magnitude  and  location  of  the  extreme  values  differ.    Under  the  HF 
approximation  the  excess  kinetic  energy  vanishes  identically  at 
absolute  zero  as  we  are  left  with  a  step  function  momentum 
distribution.    A  positive  excess  kinetic  energy  as  exhibited  by 
Pokrant  would  require  the  self-energy  beyond  the  HF  theory  (as  in  Fig. 
7.1),    That  in  turn  would  involve  a  coupling  constant  integration  of 
density  fluctuations.    This  underscores  again  the  need  for  accurate 
values  of  the  pair  distribution  function,  and  so  we  will  next 
concentrate  on  improving  the  quantal  integral  equation  results  by 
incorporating  screening  effects. 
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Table  7.1 


A  comparision  of  interaction  energies  per  particle  in  rydbergs 
for  a  density  of  rg  =  2  for  varying  temperatures  tau  =  kt/Fermi 
energy. 

-Potential /N 

Tau  Pokrant  BSS  QHNC  g(r)      QHNC  gan^r) 


.1 

.600 

.633 

.609 

.632 

.2 

.597 

.631 

.614 

.629 

.5 

.580 

.607 

.616 

.602 

.8 

.557 

.576 

.600 

.570 

1 .0 

.541 

.556 

.584 

.550 

1 .33 

.515 

.527 

.557 

.520 

1 .60 

.496 

.506 

.536 

.499 

471 
.  H  /  1 

.  t  /  O 

.  oU/ 

.4/  1 

.455 

.461 

.488 

.454 

3.20 

.412 

.417 

.440 

.410 

5.33 

.346 

.349 

.364 

.342 

8.00 

.296 

.297 

.308 

.291 

10.0 

.270 

.271 

.279 

.265 

12.5 

.245 

.246 

.253 

.241 

16.67 

.216 

.216 

.221 

.211 

20.0 

.198 

.199 

.203 

.195 

25.0 

.179 

.179 

.182 

.175 
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Table  7.2 


A  comparison  of  zero  temperature  quantum  Monte  Carlo  results  to 
those  of  Pokrant  for  the  interaction  energy  per  particle  and  excess 
kinetic  energy  for  density  values  computed  by  Pokrant. 

Ts  Quantum  Monte  Carlo  Pokrant 

■  -Pot  dKE  -Pot  dKE 


1.0  1.112  .07629  1.117  .089 


2.0  .5982  .04985  .602  .059 


3.39  .3739  .03379  .376  .041 
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Table  7.3 


A  comparison  of  excess  kinetic  energies  in  rydbergs  from  Pokrant 
and  the  BSS  methods  for  a  density  of  r^  =  2  for  varying  temperatures 
tau  =  kt/Fermi  energy.    The  ideal  fermion  kinetic  energy  per  Nkt  minus 
the  Maxwell  Boltzman  value  of  3/2  is  also  given. 

Tau  Pokrant  BSS  Ideal  Fermion 

dKE  dKE  KE/Nkt  -3/2 


.1 

.060 

.223 

4.742831e  0 

.2 

.058 

.212 

1.957581e  0 

.5 

.046 

.198 

5.432847e-l 

.8 

.031 

.103 

2.734860e-l 

1.0 

.023 

.0811 

1 .967396e-l 

1.33 

.013 

.0570 

1 .2887646-1 

1 .60 

.007 

.0428 

9.788451 e-2 

2.0 

.002 

.0286 

7.017671e-2 

2.29 

-.001 

.0220 

5.732930e-2 

3.20 

-.005 

.00967 

3.476111e-2 

5.33 

-.008 

-.00100 

1.619181e-2 

8.0 

-.008 

-.00509 

8.810031e-3 

10.0 

-.007 

-.00473 

6.305060e-3 

12.5 

-.007 

-.00523 

4.512104e-3 

16.67 

-.006 

-.00454 

2.930145e-3 

20.0 

-.005 

-.00432 

2.229816e-3 

25.0 

-.004 

-.00407 

1 .595610e-3 
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Table  7.4 


A  comparison  of  excess  kinetic. energies  (dKE  =  KE  -  KE^deal) 
and  chemical  potentials  (du  =  u  -  u^^eal)  in  rydbergs  from  a 
Hartree-Fock  decoupling  of  the  one  particle  Green's  Function  and  that 
of  Pokrant  for  a  density  of      =  0.5  at  varying  temperature,  tau  = 
kt/Fermi  energy.    Also  included  are  values  of  the  kinetic  energy  and 
chemical  potential  of  the  ideal  fermion  gas  in  rydbergs. 


Tau 

Pokrant 
dKE 

HF 
dKE 

HF 
du 

Ideal 
KE 

Ideal 
u 

.025 

-.010 

-2.442 

8.862 

14.725 

.05 

.126 

-.034 

-2.439 

8.930 

14.702 

.1 

.124 

-.112 

-2.425 

9.197 

14.610 

.2 

.103 

-.312 

-2.362 

10.188 

14.211 

.3 

-.476 

-2.246 

11.598 

13.474 

.4 

-.575 

-2.098 

13.252 

12.381 

.5 

-.020 

-.620 

-1 .941 

15.052 

10.948 

.6 

-.632 

-1.789 

16.945 

9.206 

.7 

-.624 

-1.650 

18.901 

7.184 

.8 

-.095 

-.605 

-1 .524 

20.903 

4.908 

1.0 

-.119 

-.556 

-1.313 

24.998 

-0.316 

1 .33 

-.136 

-.464 

-1 .058 

31 .917 

-10.585 

1 .60 

-.139 

-.418 

-0.909 

37.665 

-20.247 

2.00 

-.137 

-.353 

-0.750 

46.266 

-36.264 

2.29 

-.132 

-.317 

-0.664 

52.541 

-48.953 

2.67 

-.125 

-.278 

-0.577 

60.797 

-66.753 

3.20 

-.116 

-.237 

-0.489 

72.356 

-93.489 

4.00 

-.103 

-.194 

-0.395 

89.863 

-137.363 
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Figure  7.1 


Average  momentum  distribution  in  the  electron  gas  as  calculated 
from  the  random  phase  approximation,  for  =  2,  3,  4,  and  5  at  zero 
temperature. 


CHAPTER  VIII 
RESULTS  OF  THE  EXTENDED  QHNC  EQUATIONS 


In  this  chapter  we  will  consider  results  from  two  quantal 
extensions  of  classical  integral  equations.    These  were  derived  in 
Chapter  III  and  both  incorporate  the  effects  of  screening  the  long- 
ranged  coulomb  interactions.    Like  the  QHNC  method,  in  the  appropriate 
limit,  they  also  reduce  to  the  classical  HNC  equations. 

One  system,  which  we  shall  call  the  Quantal  Hartree  system,  is 
based  on  an  effective  potential  of  the  form 

-  ^^eff^^^ 

e  =  n*(r  I  VM)/n^ 

H  0 

n*(r|v^)  is  the  density  distribution  of  ideal  fermions  in  the 
presence  of  an  external  Hartree  screened  potential 

V^(r)  =  V(r)  +  p  J  d?^  V(  I  ?  -  r'  I  )  {g^"(?')  -  1} 

v(r)  is  the  bare  interparticle  potential  in  the  interacting  system 
(coulombic).    The  analogue  pair  distribution  function  n(r|v)/n 

0 

which  appears  in  the  convolution  of  the  Hartree  potential  is  to  be 
solved  self-consistently  with  the  equation 
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together  with  the  classical  Ornstein-Zernike  equation  in  terms  of 
analogue  distribution  functions.    Notice  that,  unlike  the  classical 
OCP  potential  or  the  QHNC  effective  potential,  that  which  appears  here 
is  short-ranged.    Instead  the  long-ranged  nature  of  the  coulombic 
many-body  problem  appears  in  the  analogue  bridge  function 

B*"(r)  =  F"^{PXq(0)  F{V^(r)  -  V(r)}} 

and  implicitly  in  the  OZ  relation.    A  misnomer,  the  analogue  bridge 
function  does  not  correspond  to  a  quanta!  generalization  of  classical 
bridge  diagrams.    However  it  appears  in  various  integral  equations  in 
the  same  manner  as  the  classical  bridge  function  would,  hence  its 
name.    For  example,  there  is    the  Hartree  screened  Quantal 
Percus-Yevic  equation 


g^"(r)  =  e~^^'^^^'^  {1  .  Y^"(r)  +  B^"(r)} 


We  note  in  passing  that  both  the  effective  potential  and  analogue 
bridge  function  remain  finite  at  zero  temperature. 

The  other  system  of  equations  under  consideration  we  shall  call 
the  quantal  Zwanzig  system.    It  is  based  on  an  effective  potential 

-3V^..(r) 

e  =  n*(r  I  V,)/n^ 

z  0 

wherein  the  externally  imposed  Zwanzig  potential  the  bare  interaction 

under  the  convolution  is  replaced  with  the  analogue  pair  correlation 
functionl39.162 

V^(r)  =  V(r)  +  p  ;  dr'  C^"(  I  r  -  ?')  h^"(r') 
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The  analogue  pair  distribution  function  is  solved  self-consistently 
with  the  OZ  relation  and  the  closure 

an  Xff^^^ 
g^"(r)  =  e 

This  system  of  equations  has  been  studied  extensively  at  zero 

1  98 

temperature  by  Chihara.        Because  the  effective  potential  is 
short-ranged  the  coulombic  influence  on  the  solutions  enters  only 
implicitly  through  the  OZ  equation.    The  long-range  nature  of  the 
analogue  pair  correlation  and  nodal  functions  must  be  inferred  from 
the  fact  that  the  Zwanzig  closure  can  be  written  in  the  same  form  as 
the  Hartree  closure,  except  that  the  Zwanzig  screened  potential 
appears  in  the  definition  of  the  analogue  bridge  function,  while 
ignoring  the  fact  that  the  nodal  and  bridge  functions  cancel 
identically. 

Method  of  Solution 

Both  the  Hartree  and  Zwanzig  equations  are  self-consistent  field 
calculations  for  the  distribution  of  interacting  electrons  in  the 
presence  of  an  elementary  test  charge;  both  were  solved  using  a  Picard 
iteration  by  mixing  input  and  output  to  accelerate  convergence. 
However,  there  are  subtle  differences  which  preclude  concurrently 
reviewing  their  methods  of  solution.    We  will  concentrate  on  the 
Hartree  solution  first. 

The  Hartree  system  converged  conditionally  on  the  quality  of  the 
initial  guess  for  the  analogue  pair  distribution.    Occasionally  the 
solution  at  one  temperature/density  was  an  unacceptable  seed  for 
another.    To  construct  an  initial  guess  ab  initio  we  imposed  the 
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constraint  that  the  electrons  redistribute  themselves  so  as  to 
completely  shield  the  test  charge.    It  is  found  that  many-body  effects 
screen  the  bare  coulomb  potential  into  an  exponential  damped  form^^ 

e?  -  r/d 
r  ^ 

where  the  generalized  screening  length 

(8.1) 


\  =  4Tr  e^  ^  n(vi) 
d2  dp 


reduces  to  the  Debye  screening  length  in  the  classical  limit 
^  =  (4Tr  p  e^/kT)^^^ 

and  the  Thomas-Fermi  screening  length 


1         •J    1/6     -  2 
_L  _  ,3^x         4  me 

in  the  completely  degenerate  limit.    The  ansatz  was  made  that  h*"(r) 
could  be  modeled  by  the  distribution  of  non-interacting  fermions  in 
the  Thomas-Fermi  approximation  of  a  modified  Yukawa  potential: 

a(r)  is  a  trial  function  of  the  form  l-exp(-const*r) .    This  form 
results  in  divergences  at  small  r  of  the  potential  which  mimic  the 
multiple  commutators  of  the  coulombic  1/r  pole  with  the  momentum 
operator  found  in  Wigner-Kirkwood^^  corrections  to  the  T-F 
approximation.    The  constant  was  chosen  to  ensure  total  screening, 
that  is 

ph^"(0  =  0)  =  J  dr  ph^"(r)  =  -  1 
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A  plot  of  an  initial  guess  for  h^"(r)  constructed  in  this  manner 
(for  r^  =  2  and  kt  =  1  Fermi  energy)  is  presented  as  an  illustration 
in  Fig.  8.1 . 

It  is  surprising  that  this  ansatz  is  satisfactory,  because 
asymptotically  it  is  of  the  form  (dashed  curve  of  Fig.  8.1) 

h^\r)  ~  -  -4  e- 
4irpd'^ 

(even  at  zero  temperature)  whereas  physically  h^"(r)  exhibits 
Friedel  oscillations^"^  at  zero  temperature  of  the  form 

Ze       2g       cos  (2KpX)  _  ^ 

~  ir  2         ^  ^  -       ?        9  (8.2) 

(4  +  V         x-^  2Kj-^  d^P^ 

These  arise  from  the  sharp  Fermi  surface;  it  its  not  possible  to 
construct  a  smooth  function  out  of  a  restricted  set  of  wavevectors 
q  <  k^. 

This  explanation  does  not  depend  on  the  strength  of  the  external 

perturbation,  so  the  results  should  be  implied  also  by  calculations 

based  on  the  dielectric  function  method  in  linear  response. 

There  it  is  often  claimed  that  the  oscillations  are  a  consequence  of 

the  singularity  of  the  dielectric  function  at  q  =  2k^,  but  (except 

in  the  limit  of  very  large  radius)  this  is  not  strictly  true  because 

the  slightest  blurring  of  the  Fermi  surface  due  to  thermal  excitation 

is  enough  to  destroy  the  singularity  but  does  not  remove  the 
73  200 

oscillations.    '      .    (At  large  radius  R  damping  of  the  oscillations 
requires  a  blurring  over  a  range  of  k  of  order  pi/2R-see  eq.  8.2.)  In 
fact  thermal  effects  result  in  an  exponential  damping  of  the 
oscillations. 
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That  the  starting  ansatz  remains  good  down  to  zero  temperature 
reflects  the  unimportance  of  the  asymptotic  form  in  r-space  relative 
to  a  general  agreement  in  momentum  space  with  the  actual  solution. 

The  iterative  cycle  to  solve  for  h*"(r)  is  described  as 
follows.    Given  an  initial  guess  for  h^"(r)  the  Hartree  potential  is 
computed  in  Fourier  space  as  using  the  OZ  relation 

an  ^COUL^Q^ 
1  -  p  C  (Q) 

the  coulomb  potential  is  then  expressed  as  the  sum  of  long-  and  short- 
ranged  pieces  (V^  =      +        ^^^^  possess  analytic  transforms 

p2  2 
Vg(r)  =      erfc  (ar)  Vg(r)  =  ^  erf  (ar) 

In  the  decomposition  alpha  is  an  arbitrary  parameter  to  be 
conveniently  chosen  in  a  manner  described  later.    This  allows  the 
Hartree  potential  to  be  similarly  decomposed  into  a  sum  of  the  short- 
ranged  and  "extra"  potentials  (V^  =  where 

V^(Q)  =  [pC*"(0)  V^(0)  +  Vg(0)]  /  [1  -  pC^"(Q)] 

This  extra  potential  is  transformed  numerically  into  real  space  for 
use  in  solving  the  Schroedinger  equation  in  the  Kohn-Sham  system  of 
equations. 

Fiqures  8.2  and  8.3  are  plots  of  "extra"  potentials  at  a 
temperature  of  kt  =  one  Fermi  energy  at  densities  r^  =  1.0  and 
2.00.    These  were  computed  using  the  solution  values  of  h^"(r).  It 
can  generally  be  said  that  higher  r^  values  decay  more  slowly;  at 
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the  same  temperature  all  have  comparatively  the  same  maximum  value 
which  decreases  with  higher  temperature. 

The  Hartree  potential  at  constant  density  for  varying 
temperatures  is  plotted  in  Fig.  8.4;  the  only  major  difference  being 
that  higher  r^  values  decay  more  slowly.    Plots  of  the  total  Hartree 
potential  at  varying  densities  for  the  same  temperature  are  presented 
in  Fig.  8.5.    The  interesting  feature  is  that  the  plots  are  all  of  the 
exact  same  shape—however  the  x-axis  was  scaled  by  the  temperature- 
dependent  generalized  screening  length.    Again  it  must  be  emphasized 
that  this  is  a  feature  of  the  solution  h^"(r). 

The  effective  potential  entering  into  the  system  of  integral 
equations  is  related  to  the  logarithm  of  the  density  distribution  of 
non-interacting  fermions  in  the  presence  of  the  Hartree  potential. 
Unlike  the  QHNC  effective  potential  this  has  no  analytic  solution,  and 
we  must  resort  to  a  numerical  solution  of  the  Kohn-Sham  system  of 
equations.    Computationally  this  is  a  very  expensive  process,  for 
although  n*(r|v^)  soon  reaches  its  TF  asymptotic  limit,  the 
continuum  wavefunctions  upon  which  it  is  built  must  be  calculated  to 
much  larger  radial  values  in  order  to  ensure  proper  normalization. 

An  alternative  approach  based  on  the  work  of  March  et 
JO, 203  ^5  receiving  renewed  interest, deserves 

comment.    The  expession  n*(r|v)  can  be  exactly  given  as  the  energy 
integral  over  the  Fermi  distribution  of  a  function  that  satisfies  a 
third  order  differential  equation  in  r.    Current  applications  satisfy 
the  boundary  conditions  in  terms  of  free  field  solutions  via  a 
perturbative  solution  in  orders  of  the  potential.    The  first  order 
solution  is  the  (zero  temperature)  Mott  equation^^^ 
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n 


^      (2Kp  I  r  -  s  I) 


I  r  -  s  |2 


where 


-2 


T^(x)  =  X  ^[SINx  -  xCOSx] 

This  can  be  shown  to  give  negative  values  near  the  radial  origin.  For 
such  reasons  the  (albeit  time-consuming)  Kohn-Sham  equations  were 
preferred  for  their  overall  robustness.    (For  example  the  algorithms 
for  solving  the  Schroedinger  equation  are  well  known.) 

Figures  8.6-8.10  present  solutions  of  the  KS  equations  at  zero 
temperature  for  varying  densities.    The  results  merge  with  the 
asymptotic  TF  approximation  at  smaller  radii  with  decreasing  r^ 
values,  or,  from  Figs.  8.11-8.13,  with  increasing  temperature  (as 
expected).    The  corresponding  effective  potentials  to  Figs.  8.6-8.13 
are  presented  in  Figs.  8.14  and  8.15.    As  only  the  ratio  of  the 
effective  potential  over  kt  remains  finite  at  zero  temperature,  an 
analogue  to  a  classical  potential  was  defined  by  multiplying  the  ratio 
by  1  rydberg.    We  see  that  with  decreasing  r^  the  effective 
potential  weakens  and  becomes  more  short-ranged.    The  effective 
potential  is  only  weakly  dependent  on  temperature. 

The  analogue  bridge  function,  which  arises  as  a  consequence  of 
screening,  is  also  obtained  from  the  initial  h^"(r)  via  the  Hartree 
screened  potential.    Its  calculation  is  performed  by  writing 


B°"(0)  =  3x^(00)  [V^(Q)  -  V(,Q^L.^O)] 
=  |3Xq(0.  0)  [V^(0)  -  Vg(0)] 

=  {3Xq(Q,  0)  Vj^(Q)  -  (3Xq(Q.  0)  -  3x^(00))  Vg(Q)}  -  Px(OO)  Vg(Q) 
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This  splits  the  bridge  function  into  short-  and  long-ranged  pieces. 
The  short-ranged  bridge  function  is  numerically  transformed  into  two 
parts;  the  first  part  involving  the  static  Lindhard  function  with  the 
"extra"  potential,  the  second  part  the  difference  of  the  static 
Lindhard  function  from  its  zero  wavevector  value  times  the  long-ranged 
potential.    In  the  course  of  this  decomposition  each  individual  part 
is  rearranged  into  a  machine  finite  form.    A  typical  output  for  the 
calculation  of  the  components  of  the  short-ranged  bridge  function  is 
illustrated  by  Fig.  8.16  for  r^  =  2  at  zero  temperature.  The 
corresponding  short-ranged  bridge  function  is  shown  in  Fig.  8.17. 

The  long-range  piece  of  the  bridge  function  is  identical  to  the 
long-range  tail  of  the  direct  correlation  function.    This  is  easily 
demonstrated  as  follows.    The  Fourier  transform  of  the  starting  ansatz 
for  h*"(r)  yields  for  small  wavevectors 

ph    (Q)  =  -   r 

1  +  (Qd)*^ 

which  implies  through  the  Ornstein-Zernike  relation 

pC^"(Q)  =  -  — 

(Qd)^ 

and  so  at  large  radii 

C^"(r)  ~  -  e^  Px(OO) 

The  arbitrary  constant  appearing  in  our  long-range  reference  function 
is  chosen  such  that  a  short-ranged  direct  correlation  function 

pC^"(Q)  =  p  C^"(0)  +  p3Xq(00)  Vg(0)  (8.3) 
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is,  besides  being  finite  and  short-ranged,  relatively  small  in 
magnitude  in  order  to  reduce  numerical  errors  that  occur  in  numerical 
Fourier  transformations.    We  found  that  taking  the  constant  equal  to 
the  screening  length  (see  eq.  8.1)  divided  by  the  square  root  of  two 
adequately  met  this  purpose.    As  an  illustration  the  (solution)  short- 
ranged  pair  correlation  function  for  r^  =  2  at  zero  temperature  is 
plotted  in  Fig.  8.18.    The  value  at  the  origin  corresponds  to 
approximately  one  tenth  the  value  of  the  long-range  contribution. 

With  the  effective  potential  and  analogue  bridge  functions 
specified,  the  OZ  and  closure  equations  are  then  solved  iteratively 

for  a  new  short-ranged  direct  correlation  function  in  the  manner 
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described  by  Ng.        Then  the  solution  for  the  new  (analogue)  short- 
ranged  direct  correlation  function  serves  as  a  new  input  for 
calculating  a  new  Hartree  screened  potential.    The  cycle  is  repeated 
until  the  rms  value  of  the  difference  of  input  and  output  short-range 
c^"(r) 


{  J    dr  I  c^UT.c^n  ,2jl/2 


is  smaller  than  the  rms  magnitude  of  c^"(r)  by  a  factor  of  one  part 

in  ten  to  the  fifth  power. 

The  iteration  solely  in  terms  of  short-ranged  functions  avoids 

the  divergent  results  that  plague  the  self-consistent  solution  of  the 

KS  equation,  due  to  the  long-range  nature  of  the  coulomb 
206—207  121 

interactions.  '        To  date  several  other  procedures  have  been 

proposed  to  avoid  numerical  instabilities;  for  example  Popovic 
et  al.^°^  and  Zaremba  et  al.^^^  have  used  a  parameterized 
electrostatic  potential  to  fulfill  the  Friedel  sum  rule,  while 
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Manninen  et  al .  have  made  use  of  the  Poisson  equation  substructed 
contribution  of  the  Thomas-Fermi  potential.    The  procedure  followed 

1  00 

here  is  based  on  the  work  of  Chihara. 


Ouantal  Hartree  Results 


The  solution  analogue  pair  correlation  function  obtained  from  the 
quantal  system  of  equations  was  used  to  construct  the  local  field 
correction  factor  in  the  manner  developed  in  Chapter  V.    This  was  used 
to  calculate  the  physical  pair  distribution  function  as  described  in 
Chapter  VI. 

The  results  represent  a  significant  improvement  over  the  QHNC 

data  reported  in  Chapter  VII.    Fig.  8.19  presents  the  zero  temperature 

interaction  energy  per  electron  for  six  densities  between  r  =1.0 

s 

and  3.39.  Also  plotted  are  the  interaction  energies  obtained  from  the 
high-density  expansion  of  the  correlation  energy. ^^^"^^ ^ 

^CORR  °  '^^^^  ^"  ^s  "  -^^^  ^  -^^^  ''s  ^"  ^s 
and  the  low  density  expansiontl 9] 
.876      2.65  2.94 
s  s 

The  range  of  densities  investigated  was  selected  to  bridge  the 
intermediate  density  gap  between  these  known  results. 

The  quantum  Monte  Carlo  results  do  not  differ  visually  on  the 
graph  and  a  comparison  of  Quantal  Hartree  with  MonteCarlo  values, 
along  with  those  provided  by  Pokrant,  is  presented  in  Table  8.1.  The 
quantal  Hartree  values  have  a  maximal  deviation  of  1.6  percent  at  the 
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high-  and  low-density  ends,  being  above  and  below  the  Monte  Carlo 
results,  respectively.    The  three  points  supplied  by  Pokrant  are 
uniformly  one  half  percent  above  the  Monte  Carlo  data. 

Figures  8.20-8.25  present  the  analogue  and  pair  distribution 
functions  corresponding  to  the  data  of  Table  1.    Referring  to  Fig.  26, 
we  see  that  the  pair  distribution  function  remains  positive  at  the 
radial  origin— unlike  the  RPA  or  Hubbard  theory.    A  comparison  with 
the  STLS  results  can  be  made  through  their  g(r)  values  at  the  origin: 


STLS  Hartree 


1.0 

.24 

.23 

2.0 

.04 

.18 

3.0 

.04 

.15 

The  Hartree  results  decay  smoothly  with  r^,  and  remain  much  larger 
in  magnitude.    Results  from  the  quantum  Monte  Carlo  approach  or  that 
of  Pokrant  are  not  available  for  comparison. 

The  quality  of  the  analogue/pair  distribution  functions  can  be 
expressed  quantitatively  by  the  degree  to  which  each  fulfills  the 
requirement  of  perfect  screening.    For  the  density  distribution  in  the 
presence  of  an  external  test  charge 

;  n(r  I  v)  -  n^  dr   =    p  ;  dr  h*"(r) 

=  ^^"^  ^ — V~  - 1}  =  - 1 

O-o    1  -  pC^"(0) 
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due  to  the  OZ  relation  because  the  analytic  long-range  piece  of 
c^"(q)  diverges  at  the  wavevector  origin.    Therefore  any  deviation 
from  this  value  is  a  measure  of  the  numerical  accuracy  of  results. 
The  same  cannot  be  said  of  the  physical  pair  distribution.  In 
addition  deviations  can  be  attributed  to  the  range  of  validity  of  the 
local  field  correction  factor  assumed.    Moments  for  both  h(r)  and 
h^"(r)  were  computed  in  real  space  and  are  presented  in  Table  8.2. 
Except  for  the  two  largest  r^  values  the  analogue  moment  was 
fulfilled  to  within  hundredths  of  a  percent,  the  physical  to  under  one 
percent.    This  verifies  both  the  accuracy  of  the  numerical  procedure 
and  the  validity  of  the  LFCF.    The  degradation  of  the  two  higher  r^ 
values  can  be  attributed  to  coarser  grid  sizes  precipitated  by  slowly- 
decaying  effective  potentials.    This  phenomenon  was  also  encountered 

1  QR 

and  explained  in  terms  of  phase  shift  calculations  by  Chihara. 

The  variation  of  the  interaction  energy  per  electron  with 
temperature  is  presented  in  Figs.  8.27-8.29  where  comparisons  are  made 
with  data  at  the  densities  supplied  by  Pokrant.    The  most  striking 
feature  is  the  hump  that  occurs  near  the  temperature  origin,  which 
becomes  relatively  more  pronounced  at  the  higher  r^  values. 

Because  this  feature  1)  does  not  appear  if  one  were  to  compute 
the  interaction  energy  from  h^"(r)  (see  Fig.  8.30),  and  2)  is  also 
exhibited  in  the  QHNC  equations,  its  existence  may  lie  in  the 
approximation  for  the  LFCF  assumed.    This  will  be  examined  in  the  next 
chapter. 

The  variation  of  the  pair  distribution  function  with  temperature 
is  illustrated  by  Figs.  8.31-8.34  for  a  density  of  r^  =  2.0.    At  the 
low  temperature  side  of  the  hump,  the  physical  pair  distribution 
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function  lies  above  the  analogue  distribution  (see  also  Fig.  8.22 
which  is  the  zero  temperature  plot),  at  the  peak  of  the  hump  they  are 
nearly  coincident,  while  above  the  hump  the  physical  distribution  lies 
below  the  analogue. 

Table  8.3  presents  a  comparison  of  the  analogue/physical 
screening  moments  with  temperature.    It  is  seen  that  the  analogue 
function  fulfills  the  perfect  screening  criterion  to  within  thousanths 
of  a  percent  all  the  way  down  to  zero  temperature.    The  physical 
distribution  retains  this  accuracy  at  the  high  temperature  end,  but 
suffers  a  comparatively  large  degradation  at  the  low  temperature  end 
below  the  hump.    These  two  features  reinforce  the  idea  that  the  hump 
is  an  unrealistic  feature  inherent  in  our  present  approximation  of  the 
LFCF. 

It  is  also  important  to  note  that  at  r  =  0  the  physical  pair 
distribution  function  goes  negative  at  intermediate  temperatures 
before  resuming  positive  values  with  the  advent  of  a  classical 
description.    The  onset  of  such  a  degradation  occurs  at  smaller 
temperatures  for  higher  values  of  r^.    Thus,  one  of  the  criteria 
used  to  judge  the  success  of  zero  temperature  theories,  namely  the 
positivity  of  the  pair  distribution  function,  may  be  an  illusury 
feature  with  regard  to  a  more  complete  temperature  dependent  theory. 

The  Zwanziq  Equation 

In  the  Zwanzig  system  of  equations  we  subtract  out  a  long-range 
reference  function  that  (unlike  the  quantal  Hartree  system)  does  not 
possess  an  analytic  transform  (compare  with  eq.  8.3). 
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pC^"(0)  =  pC^"(0)  +  p(JXq(0.  <!>)  Vg(Q) 

Furthermore,  instead  of  iterating  in  terms  of  c^"(r),  we  define  new 
functions  (of  dimension  rydbergs)  that  correspond  to  the  analogue  pair 
and  nodal  functions 

'^^^  -  (5Xq(Q,  0)  ^^^^  -  PXq(Q,  0) 

along  with  corresponding  definitions  involving  the  short-ranged 
analogue  decompositions.    The  iterative  solution  starts  with  an 
initial  guess  for  the  short-ranged  c(bar)  in  wavevector  space,  from 
which  we  obtain  the  short-ranged  quasi  nodal  function  by  the  OZ 
relation 

C(0)  Px^(0.  0)  C  (Q)  -  p  V  (Q) 

t  CO)  =  ^  

s^^^  1  -  3Xq(Q,  0)  C(0) 

When  the  screened  Zwanzig  potential  is  resolved  into  short-ranged  and 
excess  parts,  the  excess  potential  is  given  explicitly  in  terms  of  the 
short-ranged  quasi  nodal  function 

V^CQ)  =  V^(Q)  -  r^(Q)/p 

This  is  numerically  transformed  to  real  space  to  solve  the  Kohn-Sham 
equations.  This  results  directly  in  a  new  h^"(r),  without  having  to 
involve  a  classical  fluid  integral  equation  routine. 

The  iterative  loop  is  closed  by  taking  the  Fourier  transform  of 
the  new  h*"(r)  and  the  old  short-ranged  quasi  nodal  function  to 
obtain  a  new  guess  for  the  short-ranged  c(bar)  from  the  definition 
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The  iteration  continues  until  the  rms  difference  of  the  input  and 
output  short-ranged  c(bar)  in  real  space  is  within  an  acceptable 
tol erance. 

Like  the  quantal  Hartree  system  the  rate  of  convergence  depends 
on  the  quality  of  the  first  guess.    Unlike  the  Hartree  system  it  was 
found  that  in  lieu  of  solutions  from  other  temperature  and  density 
points  the  best  starting  input  was  a  null  function. 

One  subtle  difference  in  the  Hartree  and  Zwanzig  systems  arises 
from  the  fact  that  distribution  functions  are  of  a  necessity 
transformed  from  real  to  wavevector  space  and  back  again  in  the 
iterative  process.    This  occurs  less  frequently  in  the  Zwanzig 
system—no  classical  hnc  set  of  equations  must  be  solved  in  the 
overall  iterative  loop— and  so  one  possible  source  of  numerical  error 
is  eliminated. 

Fourier  transforms  were  accomplished  using  a  Fast  Fourier 
Transform  algorithm.    Due  to  discretization  and  wrap-around  effects 
the  asymptotic  form  of  an  FFT  transform  differs  from  an  analytic 
transform,  and  taking  inverse  FFT  transforms  of  functions  dependent  on 
analytic  transforms  may  yield  slight  errors.    Such  an  operation  occurs 
only  once  in  each  system  (involving  a  distribution  function  and  the 
analytic  Lindhard  function).    However  calculating  the  pair 
distribution  function  through  the  LFCF  requires  c(bar)  in  wavevector 
space  as  input;  this  is  supplied  by  the  Zwanzig  system  directly,  the 
Hartree  system  requires  another  mixed  transform. 
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It  is  amazing  that  in  spite  of  the  overt  difference  of  the  two 
systems  of  equations  or  the  subtle  differences  mentioned  above,  the 
results  obtained  by  independent  calculation  of  the  Zwanzig  system 
agreed  with  those  obtained  by  the  Hartree  system  within  three 
significant  figures  for  both  the  analogue  and  physical  pair 
distribution  functions  at  all  temperatures  and  densities  presented. 
This  can  only  be  ascribed  to  the  fact  that  both  systems  of  equations 
were  derived  from  the  Percus  generating  functional  method  on  the  basis 
of  a  common  quantal  bootstrap  ansatz— that  is,  the  same  nonlocal 
density  approximation  to  the  exchange  correlational  functional  in 
density  functional  theory  (refer  to  Chapter  III.)    That  is  in  the 
regions  of  temperature  and  density  considered,  the  main  effects  were 
pointedly  quantal  and  arose  in  the  same  physical  yet  mathematically 
inequivalent  ways. 

Blending 

In  classical  statistical  mechanics  the  analysis  of  the  density 
expansion  of  g(r)  has  led  to  two  "exact"  equations  which,  however, 
involve  three  unknowns: 

h(r)  =  c(r)  +  p  J  dr'  h(  I  r  -  r'  I  C(r') 
g(r)  =  exp{-  pV(r)  +  h(r)  -  C(r)  +  B(r)} 

The  first  of  these  is  the  Ornstein-Zernike  relation  defining  the 
direct  correlation  function  c(r)  in  terms  of  h(r)  =  g(r)  -  1.    In  the 
second,  a  closure  equation,  V(r)  is  the  pair  potential  and  B(r)  is  the 
sum  of  "bridge"  or  "elementary"  graphs  in  the  diagrammatic  analysis  of 
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two-point  functions.    Though  the  same  analysis  leads  to  a  formal 
relationship  between  B(r)  and  g(r),  it  involves  an  infinite  sum  of 
highly-connected  diagrams  and  so  cannot  be  utilized  in  practice.  This 
apparent  absence  of  a  simple  functional  connecting  of  B(r)  to  g(r) 
prevents  this  scheme  from  being  fully  closed. 

The  various  integral  equations  derived  in  Chapter  III  were  an 
attempt  to  derive  the  majority  of  the  effects  of  the  bridge  function 
by  starting  from  physically  sensible  generating  functionals,  the 
remainder  of  the  effects  being  lost  in  the  truncation  of  the 
generating  functional  expansion. 

In  the  practice  of  classical  statistical  mechanics  all  that  one 
can  do  to  obtain  a  more  realistic  g(r)  is  to  take  two  (approximate) 
integral  equations  that  bracket  known  (computer  simulation)  results, 
and  infer  from  their  form  what  the  true  bridge  function  must  look 
like.    The  only  guidance  is  provided  by  the  fact  that  if  we  had  the 
true  g(r)  then  the  same  equation  of  state  should  result  from  both  the 
virial  (pressure)  and  compressibility  equations.    (Based  on  scaling 
arguments  one  can  show      that  the  virial  and  internal  energy 
equations  are  equivalent  for  simple  power  law  potentials.) 

A  number  of  authors  ^'  '^'^'-^  has  proposed  imposing  this 

requirement  directly  to  the  generation  of  the  integral  equation. 

Stell^^^  has  shown  that  the  results  of  Rowlinson^^^  and  Lado^^^ 

can  be  obtained  by  a  functional  expansion  method  based  on  a  linear 

combination  of  the  PY  and  HNC  approximations.    Hutchinson  and 
102 

Conkie      studied  a  functional  of  the  form 


_  5 

W[n]  =  ^  {(n(?  I  U)  e"'^^^''^)    -  n^^} 
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where  the  parameter  S  is  varied  until  consistency  is  found.    Note  that 

S  =  1  yields  the  PY  functional,  and  by  taking  the  limit  as  S  tends  to 

zero,  the  functional  for  the  HNC  approximation  is  recovered.  A 
215 

related  approach      is  to  identify  thermodynamical ly  consistent 
truncations  of  distribution  function  hierarchies. 

Thermodynamic  consistency  is  used  in  the  Reference-HNC 

pic 

approach      which  approximates  the  B(r)  with  the  bridge  function  of 

a  short-range  (reference)  potential.    Following  Rosenfeld  and 
59  60 

Ashcroft,         who  proposed  and  extensively  documented  the  view  that 

B(r)  should  be  essentially  the  same  function  for  all  potentials,  the 

reference  potential  is  viewed  as  an  adjustable  function.    In  practice 

it  is  chosen  to  be  a  hard  sphere  system  with  variable  radius,  which  is 

172  217 

optimized  by  requiring  that  it  minimize  the  free  energy,      '  a 

condition  that  greatly  increases  the  internal  consistency.    Note  that 

the  true  bridge  function  is  a  short-ranged  function  even  for  systems 

103  218 

interacting  via  the  long-ranged  coulomb  potential.  ^'^»'-'" 

In  the  quantal  regime  the  approach  we  investigated  relied  on  two 

independently  generated  integral  equations  which  are  assumed  to 

bracket  the  true  solution  and  which  can  be  spanned  by  the  inclusion  of 

some  adjustable  function. 

For  classical  systems  with  repulsive  potentials  the  PY  and  HNC 

equations  have  this  desired  property  and  the  spanning  has  been 

2 1  9 

accomplished  with  a  mixing  of  the  two  equations.        The  requirement 


of  thermodynamic  self-consistency  has  provided  accurate  reproduction 

220 

of  Monte  Carlo  results.  Recent  work  has  considered  a  blending  of 
the  PY  and  HNC  equations  in  the  form 
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g(r)  =  e   '       {1  +  } 

where  f(r)  varies  from  zero  at  the  origin  (where  the  above  equation 
reduces  to  that  of  the  PY)  to  a  value  of  unity  at  radial  infinity 
(where  the  above  equation  becomes  the  HNC  approximation).    The  virtue 
of  this  algorithm  is  that  a  simple  form  for  the  blending  function, 
such  as 

f(r)  =  1  -  e~  "  ^  (8.4) 

suffices  and  the  adjustable  parameter  used  to  achieve  thermodynamic 
consistency  was  found  to  be  approximately  independent  of  the  coupling 
parameter.    Even  more  remarkable  was  that,  for  simple  inverse  power 
potentials,  a  simple  inverse  relationship  was  found  between  the 
adjustable  parameter  and  the  power  of  the  potential. 

Quantal  extensions  of  the  classical  HNC  and  PY  equations  have  a 
qualitatively  similar  structure.    The  classical  equations  can  be 
smoothly  obtained  from  their  quantal  extensions  by  letting  1)  x*(q) 
—a  wide  function  of  wavevector  q~be  broadened  to  a  constant  value  of 
unity,  and  2)  letting   hbar  tend  to  zero  in  the  fluctuation 
dissipation  theorem.    Because  of  this  similarity  blending  was 
attempted  with  the  hybridization  of  the  Hartree  screened  quantal  HNC 
and  PY  equations: 

-|5V"  .(r)           JY^"(r)  +  B^"(r)]  f(r)  , 
g(r)  =  e  {  W  ^  ^} 

It  was  found  that  varying  the  mixing  parameter,  alpha,  appearing 
in  eq.  8.4  from  a  value  of  1000.    (Essentially  pure  quantal  Hartree) 
down  to  1.0  resulted  in  changes  in  the  interaction  energy  of  less  than 
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one  part  in  one  thousand  for  the  density  and  temperature  regime 
presented.    This  result  is  not  unexpected  in  light  of  the  Zwanzig 
system  of  equation  results. 
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Figure  8.1 

A  typical  initial  guess  that  was  generated  ab  initio  for  rc 
1 .0  and  tau  =  1 .0 


Figure  8.2 


Plots  of  the  Hartree  screened  potential  minus  an  analytical 
short-ranged  reference  potential.    The  result  is  short-ranged  and  is 
presented  at  tau  =  1.0  for  densities  of  rg  =  1.0. 
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vh(r)  minus  vs(r)  [ryd] 


Figure  8.3 

Plots  of  the  Hartree  screened  potential  minus  an  analytical 
short-ranged  reference  potential.    The  result  is  short-ranged  and  is 
presented  at  tau  =  1.0  for  densities  of  rg  =  3.39. 
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Figure  8.4 


Plots  of  radius  times  the  screened  Hartree  potential  at  constant 
density  rg  =  1.0  and  temperature  tau  =  0,  0.4,  0.8. 
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Figure  8.5 


Plots  of  radius  times  the  screened  Hartree  potential  at  constant 
temperature  tau  =  1.0  and  rg  =  3.39,  2.0  and  1.0. 


Figure  8.6 


Plots  of  the  density  distribution  of  noninteracting  fermions  in 
the  presence  of  an  external  Hartree  potential  of  zero  temperature  for 
density      =  3.39. 
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Figure  8.7 

Plots  of  the  density  distribution  of  noninteracting  fermions  in 
the  presence  of  an  external  Hartree  potential  of  zero  temperature  for 
density      =  2.5. 
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Figure  8.8 

Plots  of  the  density  distribution  of  noninteracting  fermions  in 
the  presence  of  an  external  Hartree  potential  of  zero  temperature  for 
density  rs  =  2.0. 
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Figure  8.9 

Plots  of  the  density  distribution  of  noninteracting  fermions  in 
the  presence  of  an  external  Hartree  potential  of  zero  temperature  for 
density      =  1 .5. 
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Figure  8.10 

Plots  of  the  density  distribution  of  noninteracting  fermions  in 
the  presence  of  an  external  Hartree  potential  of  zero  temperature  for 
density      =  1 .0. 
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dens  i  ty :  n( r | v )/nO 
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Figure  8.11 

Plots  of  the  density  distribution  of  non-interacting  fermions  in 
the  presence  of  an  external  Hartree  potential  at  constant  density  rc 
=  2.0  and  temperature  tau  =  0.2. 
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Figure  8.12 

Plots  of  the  density  distribution  of  non-interacting  fermions  in 
the  presence  of  an  external  Hartree  potential  at  constant  density  re 
=  2.0  and  temperature  tau  =  0.4. 
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Figure  8.13 

Plots  of  the  density  distribution  of  non-interacting  fermions  in 
the  presence  of  an  external  Hartree  potential  at  constant  density 
=  2.0  and  temperature  tau  =  0.8. 
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Figure  8.14 

Plots  of  the  analogue  potential  employed  in  the  fluid  integral 
equations  at  zero  temperature  for  densities  rc  =  3.39,  2.5,  2.0, 
1.5,  1.0,  respectively.    Only  the  potential  divided  by  kt  remains 
finite  at  zero  temperature:    For  comparison  with  classical  potentials 
we  multiplied  by  1  ryd. 
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Figure  8.15 

Plots  of  the  analogue  potential  employed  in  the  fluid  integral 
equations  at  constant  density      =  2.0  for  temperatures  tau  =  0.2, 
0.4,  0.8,  respectively.    Only  the  potential  divided  by  kt  remains 
finite  at  zero  temperature.    For  comparison  this  potential  is  that 
quantity  multiplied  by  1  ryd. 
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Figure  8.16 

The  construction  of  the  short-ranged  analogue  bridge  Function 
from  the  inverse  Fourier  transform  of  (top)  the  static  Lindhard 
function  times  the  Fourier  transform  of  the  excess  potential  and 
(bottom)  the  inverse  Fourier  transform  of  the  difference  of  the  static 
Lindhard  function  minus  its  wavevector  origin  times  the  reference 
analytic  long-range  potential.    Illustrated  for      =  2.0  and  zero 
temperature. 


Figure  8.17 


The  total  short-ranged  analogue  bridge  function  constructed  from 
the  parts  of  Fig.  8.16.  at      =  2.0  and  zero  temperature. 


Figure  8.18 


Illustration  of  a  typical  short-ranged  analogue  direct 
correlation  function  for      =2.0  and  zero  temperature. 
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Figure  8.19 

Plot  of  the  zero  temperature  results  for  the  interaction  energy 
per  particle  obtained  from  the  quantal  Hartree  equation  versus 
density.    Analytic  high-  and  low-density  asymptotic  forms  are  also 
included. 
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Figure  8.20 

Plot  of  the  physical  and  analogue  pair  distribution  function  at 
zero  temperature  for  densities  rj  =  1.0.    The  physical  pdf  was 
computed  under  a  static  LFCF. 
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Figure  8.21 

Plot  of  the  physical  and  analogue  pair  distribution  function  at 
zero  temperature  for  densities  rg  =  1.5.    The  physical  pdf  was 
computed  under  a  static  LFCF. 
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Figure  8.22 

Plot  of  the  physical  and  analogue  pair  distribution  function  at 
zero  temperature  for  densities      =  2.0.    The  physical  pdf  was 
computed  under  a  static  LFCF, 
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Figure  8.23 

Plot  of  the  physical  and  analogue  pair  distribution  function  at 
zero  temperature  for  densities      =  2.5.    The  physical  pdf  was 
computed  under  a  static  LFCF. 
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Figure  8.24 

Plot  of  the  physical  and  analogue  pair  distribution  function  at 
zero  temperature  for  densities      =  3.0.    The  physical  pdf  was 
computed  under  a  static  LFCF. 
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Flgure  8.25 

Plot  of  the  physical  and  analogue  pair  distribution  function  at 
zero  temperature  for  densities  rg  =  3.39.    The  physical  pdf  was 
computed  under  a  static  LFCF. 


-177- 


Figure  8.26 

Pair  correlation  function  in  the  degenerate  electron  plasma 
two  values  of  rg  in  the  range  of  metallic  densities  according  to 
various  theories  of  correlations. 


-178- 


Figure  8.27 


Plots  comparing  the  interaction  energies  versus  temperature 
obtained  from  the  Quanta!  Hartree  equations  along  with  those  of 
Pokrant,  at  densities  rg  =  1.0. 
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Figure  8.28 

Plots  comparing  the  interaction  energies  versus  temperature 
obtained  from  the  Quanta!  Hartree  equations  along  with  those  of 
Pokrant,  at  densities      =  2.0. 
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Figure  8.29 


Plots  comparing  the  interaction  energies  versus  temperature 
obtained  from  the  Quantal  Hartree  equations  along  with  those  of 
Pokrant,  at  densities      =  3.39. 
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Table  8.1 


A  comparison  of  interaction  energies  per  particle  at  zero  temperature 
for  varying  densities. 


Monte-Carlo 

-Potenti  al 
Hartree 

[ryd] 
%diff 

Pokrant 

%dif1 

1 .0 

1 .112 

1 .1293 

1.5 

1.117 

.56 

1.5 

.7731 

.78274 

1 .2 

2.0 

.5982 

.60021 

.34 

.602 

.64 

2.5 

.4905 

.49028 

.04 

3.0 

.4171 

.41228 

1.16 

3.39 

.3739 

.36783 

1 .63 

.376 

.56 
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Figure  8.30 

Plot  comparing  the  interaction  energy  versus  temperature  as 
obtained  from  the  pdf  along  with  that  which  would  be  obtained  from 
analogue  pdf  at  density      =  2.0.    The  results  of  Pokrant  are  also 
supplied. 
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Figure  8.31 

Plots  of  the  physical  and  analogue  pair  distribution  functions  at 
constant  density      =  2.0  for  temperature  tau  =  0.2.    The  physical 
pdf  was  computed  under  a  static  LFCF. 


-184- 


Figure  8.32 

Plots  of  the  physical  and  analogue  pair  distribution  functions  at 
constant  density      =  2.0  for  temperature  tau  -  0.4.    The  physical 
pdf  was  computed  under  a  static  LFCF. 
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Figure  8.33 

Plots  of  the  physical  and  analogue  pair  distribution  functions  at 
constant  density  rg  =  2.0  for  temperature  tau  =  0.6.    The  physical 
pdf  was  computed  under  a  static  LFCF. 
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Figure  8.34 


Plots  of  the  physical  and  analogue  pair  distribution  functions  at 
constant  density      =  2.0  for  temperature  tau  =  0.8.    The  physical 
pdf  was  computed  under  a  static  LFCF. 
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Table  8.2 


Total  screening  moments  of  the  analogue  and  physical  pair  distribution 
functions  obtained  from  the  Quantal  Hartree  equations  at  zero 
temperature  and  varying  densities. 

r^  Analogue  Physical 


1.0  -1.000207  -0.994870 

1.5  -1.000313  -1.009535 

2.0  -1.000701  -0.985722 

2.5  -0.999976  -0.992360 

3.0  -0.994812  -0.966804 

3.39  -0.992707  -0.964524 
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Table  8.3 


Total  screening  moments  of  the  analogue  and  physical  pair  distribution 
functions  obtained  from  the  quantal  Hartree  equations  at  a  density  of 
rs  =  2.0  at  varying  temperatures  tau  =  kt/Fermi  energy. 


Tail 

1  au 

Ana  1 ogue 

Physi cal 

1.0 

-0.9999538 

-0.9999548 

0.8 

-1 .000045 

-1.000018 

0.5 

-1 .000118 

-1.000069 

0.4 

-1 .0000151 

-1 .000097 

0.2 

-0.9999595 

-0.9982337 

0.0 

-1 .000699 

-0.9857216 

CHAPTER  IX 

DYNAMICS  OF  THE  LOCAL  FIELD  CORRECTION  FACTOR 


We  saw  in  the  previous  chapter  that  the  zero  temperature  results 
obtained  from  the  Quanta!  Hartree  system  of  equations  were  in  fair 
agreement  with  the  Quantum  Monte  Carlo  results,  but  at  finite 
temperatures  an  anomalous  hump  in  the  interaction  energies  per  particle 
was  observed  near  the  temperature  origin.    This  feature  was  found  to  be 
exhibited  by  other  systems  of  integral  equations.    It  was  not  present 
if  we  calculated  the  interaction  energies  with  the  analogue  pair 
distribution,  and  was  accompanied  by  a  degradation  in  the  value  of  the 
perfect  screening  moment.    These  features  indicate  a  breakdown  of  the 
validity  of  the  static  local  field  correction  approximation  made  in 
Chapter  V.    In  the  present  chapter  we  will  investigate  a  refinement  of 
the  LFCF  which  incorporates  dynamic  effects. 

As  our  starting  point  we  will  again  consider  the  auto  correlation 
of  density  fluctuations 

C(Q,  t)  =  ^  (n^Ct)  I  nq(o)) 

That  is  the  "Kubo  function"  in  the  scalar  product  notation  presented 
in  Chapter  V.    Recalling  the  definition  of  the  Liouvillian  time 
evolution  operator  and  the  identity 
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(LFIG)  =  -  I  <^  [F,  G*]>  (9.1) 

we  find  that  the  Kubo  function  is  related  in  real  time  to  the  density 
response  function  as 

^  C(Q,  t)  =  -  zi  x"(0.  t) 
while  the  dynamic  and  static  physical  responses  are  respectively  given  by 
x(Q.  t)  =  e(t)  ^  C(Q.  t)  (9.2) 

x(0)  =  x(0.  z  =  (0  +  ie)  =  C(0,  t  =  0) 

We  employ  the  convention  that  tilde  denote  temporal  Laplace 
transforms.    The  response  is  defined  as  in  Chapter  I  to  be  consistent 
with  the  literature  and  not  in  the  dimensionless  form  of  Chapter  III. 

The  dynamics  of  the  Kubo  function  will  be  explored  through  the 
continuity  equation 

d^  2 

^  C(Qt)  =  -  Q'^  J(0.  t) 

dt"^ 

using  the  memory  function  formalism^^^  to  describe  the 
autocorrelation  of  the  longitudinal  current 

J^Qt)  =  ^  (Lnq(t)  I  Lnq(o)) 
As  an  autocorrelation  function  J(0,t)  satisfies  the  Mori  equation 

^  J(Qt)  +  J   df  M(0.  t  -  f)  J(0.  f)  =  0 

0 

where  the  memory  function  M(0,t)  remains  unspecified  and  t>0.  However 
its  zero  time  behavior  can  be  obtained  by  differentiating  the  Mori 
equation  once 
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h2  t  . 

^  J(Ot)  +  M(0,  t  =  0)  J(Q.  t)  +  J    dt'  J(Ot')  Jt  M(Q,  t  -  t')  =  0 

dr  0  °^ 


from  which 

E  M(Q.  t  =  0)  =  -  ^  J(Ot) 
°  df^ 


/  J(Q.  t  =  0) 

t=o 


The  large  time  behavior  of  the  memory  function  is  ascertained  by 
splitting  it  up  into  its  asymptotic  limit  plus  a  piece  that  vanishes 
at  large  times 

M(t)  =       +  M^(t) 

We  then  take  the  Laplace  transform  of  the  Mori  equation 

-  J(t  =  0)  -  iZ  J(Z)  -      4^-  +  J(Z)  M^(Z)  =  0 

(for  imaginary  z>0)  and  consider  expanding  this  equation  in  powers 
of  z;  for  example  letting 

3(1)  =  J(Z  =  0)  +  Z  J^^\z)  +  Z^  J^^^Z)  +  ... 

The  first  two  terms  can  be  specified  by  considering  the  solution  to 
the  continuity  equation  expressed  in  the  form 

2  * 

C(Qt)  =  x(Q)  -      ;    dt'  J(Qt')  {t  -  t'}  (9.3) 

0 

A  sufficient  condition  that  the  Kubo  function  vanish  at  infinite  times 
is  provided  by  the  assumption 

00  00 

J   dt  J(Qt)  =0  J   dt  t  J(Qt)  =  x(Q)/Q^ 

0  0 

Therefore,  because 
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J(Z)  M^(Z) 

must  be  of  least  order  z  (only  if  M^(q,z)  were  a  constant  would  it 
have  a  1/z  behavior  near  the  origin),  we  see  that  the  zeroth  order 
term  of  the  Laplace  transformed  Mori  equation  yields 

=  -  i  J(t  =  0)  /  J^^^Z  =  0) 

Furthermore  the  autocorrelation  of  the  longitudinal  current  can  be 
evaluated  exactly  at  t  =  0  through  eq.  9.1.    It  is  essentially  a 
manifestation  of  the  f-sum  rule: 

0(Q.  t  =  0)  =  S 
m 

We  are  now  in  a  position  to  discuss  approximations  for  the 
complex  response  function;  this  function  inserted  into  the 
fluctuation-dissipation  theorem  results  in  the  dynamic  structure 
factor.    Now  it  is  clear  from  eqs.  9.2  and  9.3  that  the  complex 
response  is  directly  related  to  the  Laplace  transform  of  the 
longitudinal  current  autocorrelation  as 

which  by  virture  of  the  Laplace  transform  of  the  Mori  equation  becomes 
x(Qz)  =  —o  2  ::   (9.4) 

"o    -  2    +         -  "o^  ^"  M2(0.z)) 

Here  we  have  introduced  a  unit  normalized  memory  function  M2(q,t)  by 
the  definition 

M(Qt)  =  M„  +  (Mq  -  M^)  M2(0t) 
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where  the  scale  frequencies  are  given  by 
"o^  =      =  «to^»/x(0) 

"oo^  E       =  (L^       I       nq)/(L       I  L  n^)  =  «to^»/«£o^» 

In  the  above  equations  we  have  employed  the  double  brackets  to  denote 
frequency  moments  of  the  density  response. 

We  recall  that  by  definition  M2(q,  t  =  0)  =  1  and  this  implies 

that 

M2(Q.  z  -  <»)  =  i 

so  that  eq.  9.4  gives  the  proper  value  of  the  fourth  frequency  moment 
of  the  density  response.    The  condition  M2(q,  t  =  infinity)  =  0 
implies  that 

-  iz  M2(Qz) 

vanishes  as  z  tends  to  zero,  yielding  the  correct  static 
compressibility.    For  intermediate  times  we  have  no  real  knowledge  of 
M2(q,t)  and  it  is  here  that  one  has  to  make  a  reasonable  ansatz.  We 
make  the  approximation  that  is  is  given  by  M2^(q,t),  the 
noninteracting  fermion  memory  function. 

A  dynamic  local  field  correction  factor  can  be  obtained  by 
inverting  eq.  9.4 

1  1   "co  1 

X(Qz)      ^^^^      «a3S>  ^^^^  2F 

and  then  subtracting  the  corresponding  equation  for  a  system  of  ideal 
fermions 
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x(Oz)     x^(Qz)  "  ^(Q)     ^0^^^      '  «co2»   "   x(Q)  -  x,(0) 

(-iz  M2p(0z)) 
In  the  above  derivation  we  have  used  the  fact  that 
«M^»  =  =  pQ^/m 

The  remaining  scale  frequencies  are  given  explicitly  by^^ 


where      is  the  ideal  energy-wavevector  relation,  Wp  is  the  plasma 

87  221 
frequency  and  • 

2 

(0^(1  -  P(Q))  =^S  dr  g(r)  [1  -  cos  Qz]  ^  V(r) 

for  neutral  fluids  with  general  non-transformable  potentials,  while 

=  co^  +  ^  ;  dr  [g(r)  -  1]  {1  -  cos  Qz}  ^  V(?) 
^  8z 

for  the  OCP  or  long  ranged  potentials. 

Several  comments  can  now  be  made  about  the  dynamic  local  field 
correction  factor  V(k,z)^^^'^^^  appearing  in  the  following  equation 


x(Kz)  =  ° 


x.(KZ) 


1  -  V(K2)  Xq(Kz) 


First,  V(k,z)  can  be  written  in  the  form 

V(Kz)  =  V^^CK)  +  [Vp^(K)  -  V3^(K)]  [-  iz  M2p(Kz)] 
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where  V^^(k)  is  the  static  local  field  correction  factor  we  have 
previously  employed,  equivalent  to  the  quantum  mechanical  direct 
correlation  function.    Secondly,  because  V^^Ck)  and  V^y(k)  are 
purely  real,  and  from  the  known  analytic  properties  of  the  dynamic 
ideal  fermion  response  (temperature  Lindhard  function),  we  see  that 


X(K,  Z  =  CO  ±  ie)  =  x'(Kto)  ±  ix"(Ka)) 

and  so  the  analysis  carried  out  for  the  static  LFCF  in  Chapter  VI 
remains  in  force,  and  we  can  calculate  the  fluctuation  dissipation 
theorem  from  values  along  the  purely  imaginary  frequency  axis. 
Thirdly,  the  dynamic  LFCF  merely  changes  the  static  form  of  V^^Ck) 
over  to  that  of  V^y^t^)  at  larger  (imaginary)  frequencies.  The 
general  dynamic  dependence  can  be  straightforwardly  derived  by  making 
the  gaussian  approximation  for  the  memory  function 

-  1  a^  t^ 
M2(Kt)  =  e  ^ 

where  in  general  "a"  is  a  function  of  wavevector  k.  (Note  that  this 
satisfies  the  necessary  temporal  boundary  conditions.)  We  find  that 
for  purely  imaginary  frequencies  z  =  i  nu  (nu  real) 

U 

-  iz  M^iKz)  =  ^  ^  I  e   ^  erfc(— ^) 


V  2a 


~  (1  -  -"2)      as      V  -  <» 

V 


=  a  V  2  as  V  =  0 

Finally  we  note  that  the  dynamically  mixed  LFCF 
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4  e^i  (<T>  -  <T>r)  5 

depends  on  the  excess  kinetic  energy  and  the  pair  distribution 
function  itself  (contained  in  P(q)).    As  this  latter  quantity  is  what 
we  hope  to  obtain  from  the  calculation,  we  see  that  this  algorithm 
must  in  principle  be  calculated  self-consistently.    In  practice  P(q) 
depends  only  weakly  on  the  pair  distribution  function  itself  and  so  we 
may  reasonably  approximate  it  by  h^"(r).    This  is  easily 
demonstrated  as  follows.    For  the  OCR  the  volume  integral  for  P(q)  can 
be  reduced  to  the  one-dimensional  radial  integral 

P(0)  =  2  ;  dr[g(r)  -  1]  J  .  - 

^  (Qr)^  (Qr)-^ 

We  note  that  the  terms  in  the  braces  are  related  to  the  Bessel 
function  of  order  5/2.    This  is  an  oscillatory  function  which  is 
essentially  damped  out  by  its  first  node,  which  occurs  at  qr  =  5.763. 
For  q  much  smaller  than  the  scale  of  h(r)  we  can  series  expand  the 
braces  and  obtain 

where  Pot  is  the  interaction  energy  per  particle 

1  J 
Pot  =  ^  p  J  dr  [g(r)  -  1]  ^ 


while  for  large  q  the  only  contribution  to  the  integral  occurs  from 
around  the  radial  origin,  and  we  may  pull  the  slow  dependence  of  h(r) 
out  from  under  the  integral  by  series  expanding  it  about  its  value  at 
the  origin.    To  lowest  order  we  obtain 
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P(0)  =  -  f  If  =  0)  -  2Mn 


r=o 


This  asymptotic  value  is  quickly  reached  and  P(q)  does  not  depend 

critically  on  the  shape  of  h(r)  but  rather  on  the  single  numbers 

h(r  =  0)  and  Pot,  which  are  fairly  well  approximated  by  h*"(r)  values. 

As  for  the  dependance  on  the  excess  kinetic  energy,  it  can  be 
shown^^  that  the  successful  (zero  temperature)  STLS  theory  can  be 
derived  in  part  by  assuming  that  its  effect  is  negligible.  This 
approximation  has  also  been  made  by  other  authors. ^^^"^^^  However, 
the  fact  that  it  would  be  the  dominant  contribution  to  V^^Cq)  at 
large  wavevector  makes  this  approximation  suspect,  and  we  shall 
discuss  the  effect  of  its  inclusion  later. 


Calculating  the  Memory  Function 


Defining  the  Lindhard  function  as  in  Chapter  III,  that  is 


Px^COco)  =  ?  ;  ^  -  ^>  -    Im  CO  >  0 


we  saw  from  the  previous  section  that  for  values  of  w"  along  the 
imaginary  frequency  axis 

{w"  M(Oco")} 

In  units  of  a  frequency  parameter  nu  defined  as 
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the  dimensionless  function  R(q,co")  =  a)"M(q,to")  can  be  calculated  as 

r       1  /fl^v    (^_±J^\\       ff^Li\  __]_\ 

^3x^(0v)     ^2m^       ,^2           ^4ni^  2    "  3x(0)^ 
R(Qv)  =   2   2~2  

{2  <T>p}  ^  -  ^} 

The  braces  group  the  inverse  of  the  static  and  dynamic  Lindhard 
function  with  the  inverse  of  their  respective  leading  asymptotic 
terms.    The  function 

PXq     2m  2 

for  a  density  of  r^  =  2.0  and  at  zero  temperature  is  presented  in 

224 

Fig,  9.1.  It  is  non-analytic  about  the  Fermi  wavevector,  at  the 
wavevector  origin  it  assumes  the  value 

f,^  «K°» 


2'"  «K-2» 

while  at  large  wavevectors  it  approaches  a  constant  value  of 

2  «K^» 

3  «K°» 

(Here  double  brackets  denote  the  spherical  wavevector  average  over  the 
Fermi  distribution).    These  limiting  forms  are  equally  applicable  at 
finite  temperatures  where  Fig.  9.2  (r^  =  2  and  tau  =  kt/Fermi  energy 
=  0.2)  illustrates  that  a  slight  increase  in  temperature  modifies  the 
cusp. 

The  braces  involving  the  dynamic  Lindhard  function  exhibit 
similar  temperature-sensitive  features  which  make  calculating  R(q,w) 
non-trivial.    Specifically,  the  small 
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V  =  40^ 
Q    +  V 

asymptotic  expansion  of  the  dynamic  Lindhard  function  quickly  diverges 
at  finite  temperatures,  furthermore,  after  cancelling  the  reciprocal 
of  the  leading  term  what  is  left  in  the  numerator  of  R(q,w)  closely 
approximates  the  denominator.    The  physically  meaningful  quantity 
R(q,w)  -  1  has  thereby  suffered  a  quick  loss  of  significant  figures. 
(It  must  be  remembered  that  the  dynamic  Lindhard  function  is 
calculated  by  a  numerical  quadrature  and  starts  out  with  only  a  finite 
number  of  significant  figures.)    The  quality  of  the  numerical  results 
for  the  memory  function  becomes  increasingly  questionable  with 
temperature. 

At  zero  temperature  the  memory  function  is  analytically 
calculable,  and  the  salient  features  are  illustrated  in  Figs. 
9.2-9.3.    The  effects  of  small  but  finite  temperatures  (kt/Fermi 
energy  =  0.2)  corresponding  to  Figs.  9.2-9.3  are  presented  in  Figs. 
9.4-9.5. 

Figure  9.2  plots  the  dimensionl ess  memory  function  (at      =  2 
and  zero  temperature)  versus  wavevector  in  bohr  for  increasing 
imaginary  frequencies  of  1 ,  2  and  3  ryd.    Higher  frequencies  exhibit  a 
higher  shoulder  at  the  Fermi  wavevector.    Finite  temperatures 
eliminate  the  shoulder  (Fig.  9.4),  and  higher  frequencies  take  on  a 
gaussian  apperance. 

The  zero  temperature  frequency  dependence  is  illustrated  by  Fig. 
9.3  for  increasing  values  of  wavevector.    Finite  temperature  effects 
(Fig.  9.5)  slow  the  convergence  to  the  asymptotic  value  of  unity  but 
otherwise  alter  no  general  features. 
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Results 


The  quantal  Hartree  system  of  equations  was  re-run  at  zero 
temperature  for  varying  densities  employing  the  dynamic  local  field 
correction  factor  while  neglecting  the  excess  kinetic  energy 
contribution.    A  comparison  of  the  obtained  interaction  energy  per 
particle  is  presented  with  the  Quantum  Monte  Carlo  results  in  Table 
9.2.    It  is  seen  that  the  results  from  the  dynamic  LFCF  are  uniformly 
one  to  two  percent  greater  than  the  Monte  Carlo  data.    This  represents 
a  degradation  in  the  quality  of  results  vis  a  vis  that  obtained  from 
the  static  local  field  correction  factor.    On  the  other  hand,  the 
quality  of  the  pair  distribution  function  as  measured  by  the  criterion 
of  perfect  screening  is  seen  in  Table  9.2  to  be  uniformly  improved. 
The  pair  distribution  functions  themselves  are  presented  in  Figs.  9.6 
and  9.10.    Overall  the  values  of  the  pair  distribution  function  are 
smaller  at  the  origin  than  they  were  under  the  static  LFCF.  These 
values  are  compared  along  with  the  analogue  values  at  the  origin  in 
Table  9.3.    Curiously  enough  the  static  LFCF  values  are  uniformly 
larger  than  the  analogue  values.    This  was  the  case  for  temperatures 
on  the  low  side  of  the  temperature  "hump"  discussed  in  Chapter  VIII. 
Under  the  dynamic  LFCF  no  such  general  statement  can  be  made. 

Finite  temperature  effects  on  the  interaction  energies  per 
particle  were  next  investigated  for  those  low  temperatures  where  the 
memory  function  could  still  be  considered  accurate.    It  was  found  that 
the  values  generally  followed  the  hump  obtained  from  the  static  LFCF 
but  with  values  shifted  by  an  amount  equal  to  that  at  zero  temperature. 
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The  inclusion  of  the  excess  kinetic  energy  contribution  to  the 
dynamic  LFCF  produced  pair  distribution  functions  that  had  undesirable 
features  near  the  radial  origin.    This  is  illustrated  in  Figs. 
9.11-9.12  at  zero  temperature  for  two  densities  supplied  with  an 
artificially  small  value  of  0.01  ryd  for  the  excess  kinetic  energy. 
More  realistic  values  based  on  quantum  Monte  Carlo  data  resulted  in  a 
greater  divergence  at  the  origin.    From  the  data  in  Table  9.4  it  is 
seen  that,  in  addition  to  negative  values  of  the  PDF  at  the  origin, 
the  interaction  energies  per  particle  suffered  in  comparison  with  the 
Monte  Carlo  data  with  no  improvement  of  the  total  screening  moment. 

These  unsatisfactory  corrections  comprise  the  main  result  of  this 

chapter:    namely,  that  a  higher  order  approximation  in  the  Mori  memory 

function  formalism,  in  such  a  way  that  the  first  and  third  frequency 

moment  sum  rule  is  simultaneously  satisfied,  is  in  fact  not  a  fruitful 

course  to  follow.    Rather,  diagrammatic  analysis  at  zero 
1 9 

temperature     has  shown  the  importance  of  modeling  the  finite 

lifetime  of  the  excited  electron  and  hole  states.    These  lifetime 

effects  are  not  taken  into  account  when  one  chooses  a  real  local  field 

correction  factor  as  was  done  here.    On  the  other  hand  past  attempts 

at  incorporating  such  lifetime  effects  within  the  framework  of  the 
14-15 

Mon  theory         have  led  to  violations  of  the  continuity  equation. 

One  should  also  note  that  the  Mori  memory  function  formalism, 
which  is  essentially  a  truncation  of  a  continued  fraction  expansion, 
can  be  viewed  as  having  a  basis  in  the  belief  that  the  high  frequency 
expansion  of  the  complex  response  might  be  possible  up  to  infinite 
order.    The  following  fact,  however,  is  worthy  of  particular 
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attention.    An  asymptotic  form  of  Imag  xCq.z)  for  a  degenerate 
electron  liquid  is  given  by  Click  and  Long^^^  as 

p  6  2 

Imag  x(Ql)  =  vCQ)     ff^  ^  q  >      „  »  E  (9.5) 

ir  a^  (mco) 

It  arises  from  two-pair  excitations  of  second-order  perturbation. 
When  substituted  into  the  fifth  moment  integral  the  result  is  quite 
divergent 

r*   .  5-11/2 

J         dtd  U     CO  -»  00 

(l) 

c 

(w^  is  the  minimum  value  of  w  for  which  eq.  9.5  holds  with 
sufficient  accuracy.)    This  divergence  would  appear  when  truncating 
(in  a  Mori  continued  fraction  sense)  at  the  next  level  beyond  that 
employed  in  this  chapter  to  derive  the  dynamic  LFCF.    A  priori  this  is 
not  sufficient  to  discount  the  accuracy  of  the  dynamic  LFCF  derived  in 
this  chapter.    With  hindsight  one  may  infer  that  the  divergence 
heralds  the  degradation  of  the  dynamic  LFCF  results  over  those 
obtained  using  a  static  LFCF  obtained  from  the  Mori  method. 
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Table  9.1 


Interaction  energies  per  electron  [ryd]  obtained  from  a  pair 
distribution  function  calculated  under  the  dynamic  Local  Field 
Correction  Factor  approximation  from  the  Quantal  Hartree  equations. 
Results  are  for  zero  temperature  with  neglect  of  excess  kinetic  energy 
contributions  to  the  LFCF.    Comparison  is  made  with  Quantum  Monte 
Carlo  results. 

-Pot 

Monte  Carlo      Quanta!  Hartree  %diff 


1.0 

1.112 

1.1358 

2.1 

1 .5 

.7731 

.79013 

2.2 

2.0 

.5982 

.60870 

1 .7 

2.5 

.4905 

.49904 

1 .7 

3.0 

.4171 

.42154 

1.0 

I 
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Table  9.2 


Total  screening  moments  for  the  pair  distribution  function  calculated 
under  the  dynamic  local  field  correction  factor  approximation  from  the 
Quantal  Hartree  equations.    Results  are  for  zero  temperature  with 
neglect  of  the  excess  kinetic  energy  contribution  to  the  LFCF. 
Comparisons  is  made  to  the  total  screening  moments  obtained  under  the 
static  LFCF. 


Dynamic 

Static 

1 .0 

-0.9955281 

-0.994870 

1.5 

-1 .007725 

-1 .00955 

2.0 

-0.9878318 

-0.985772 

2.5 

-0.9933104 

-0.992360 

3.0 

-0.9705490 

-0.966804 
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Table  9.3 


Values  of  (minus  one  times)  h(r)  =  g(r)  -  1  at  the  radial  origin  at 
zero  temperature.    Comparison  is  made  with  results  obtained  from  the 
Quantal  Hartree  equations  under  the  static  and  dynamic  local  field 
correction  factor  approximations.    The  values  corresponding  to  the 
analogue  pair  distribution  are  also  presented. 


Dynami  c 

Analogue 

Static 

1.0 

.81565 

.80388 

.77287 

1.5 

.87668 

.89783 

.80273 

2.0 

.92277 

.94465 

.81942 

2.5 

.96570 

.96961 

.82914 

3.0 

1 .01480 

.98359 

.84664 
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Table  9.4 


Results  for  the  dynamic  local  field  correction  factor  with  the 
inclusion  of  an  excess  kinetic  energy  value  of  0.01  ryd  at  zero 
temperature. 

-Pot  -screening  -h(r  =  0) 


1.0  1.1374 

1.5  .79253 

2.0  .61189 

2.5  .50291 


.9955540  .88120 

1.007454  1.06873 

.9881874  _  1.33991 

.9929164  2.38550 
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Figure  9.1 


Difference  in  the  reciprocals  of  the  static  temperature  Lindhard 
function  and  its  leading  asymptotic  form. 
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lindhard  memory  function 
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Figure  9.2 


Lindhard  memory  function  versus  wave  vector  at  zero  temperature 
varying  imaginary  frequencies. 
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Figure  9.3 


Lindhard  memory  function  versus  imaginary  frequency  at  zero 
temperature  for  varying  wave  vectors. 
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lindhard  memory  function 


[1/bohr] 


Figure  9.4 

Lindhard  memory  function  versus  wave  vector  at  finite  temperature 
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Figure  9.5 


Lindhard  memory  function  versus  imaginary  frequency  at  finite 
temperature  tau  =  0.2  for  varying  wave  vectors. 
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physical/analogue  h(r) 


Figure  9.6 

Plot  of  the  zero  temperature  physical /analogue  pair  distribution 
functions  for  rg  =  1.0.    The  physical  PDF  was  computed  under  a 
dynamic  LFCF  approximation  neglecting  excess  kinetic  energy 
contributions. 
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physical/analogue  h(r) 
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Figure  9.7 

Plot  of  the  zero  temperature  physical /analogue  pair  distribution 
functions  for      =  1.5.    The  physical  PDF  was  computed  under  a 
dynamic  LFCF  approximation  neglecting  excess  kinetic  energy 
contributions. 
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Figure  9.8 

Plot  of  the  zero  temperature  physi cal /analogue  pair  distribution 
functions  for  rg  =  2.0.    The  physical  PDF  was  computed  under  a 
dynamic  LFCF  approximation  neglecting  excess  kinetic  energy 
contributions. 
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Figure  9.9 

Plot  of  the  zero  temperature  physi cal /analogue  pair  distribution 
functions  for      =  2.5.    The  physical  PDF  was  computed  under  a 
dynamic  LFCF  approximation  neglecting  excess  kinetic  energy 
contributions. 
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Figure  9.10 

Plot  of  the  zero  temperature  physi cal /analogue  pair  distribution 
functions  for      =  3.0.    The  physical  PDF  was  computed  under  a 
dynamic  LFCF  approximation  neglecting  excess  kinetic  energy 
contributions. 
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Figure  9.11 

Plot  of  the  zero  temperature  physi cal /analogue  pair  distribution 
functions  for      =  2.0.    The  physical  PDF  was  computed  under  a 
dynamic  LFCF  approximation  including  an  excess  kinetic  energy 
contribution  of  0.01  ryd. 
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Figure  9.12 

Plot  of  the  zero  temperature  physi cal /analogue  pair  distribution 
functions  for      =  1.5.    The  physical  PDF  was  computed  under  a 
dynamic  LFCF  approximation  including  an  excess  kinetic  energy 
contribution  of  0.01  ryd. 


CHAPTER  X 
CONCLUSIONS 


In  this  dissertation  we  have  examined  the  non-zero  temperature 
quantal  electron  gas.    Except  for  the  approach  espoused  by  Pokrant^^ 
and  his  antecedents,  few  numerical  results  of  thermodynamic  properties 
at  intermediate  degeneracy  have  been  previously  published,  despite  the 
overwhelming  amount  of  work  that  has  been  conducted  on  the  fully 
degenerate  or  completely  classical  electron  gas  systems.    Of  late  this 

has  been  changing:  temperature-dependent  effects  have  been  published 

20  7fi 
by  Dandrea  et  al .     and  Tanaka  et  al .     using  a  dielectric 

formalism  of  the  many-body  problem  by  employing  the  non-zero 

temperature  Lindhard  ideal  fermion  density  response  functions. 

However  the  effects  of  electron-electron  interactions  have  been 

introduced  through  static  local  field  corrections  that 


(1)  have  employed  the  STLS  approximation  which  is  known  to  give 
worse  results  than  the  classical  HNC  equation  in  the 
non-degenerate  regime, 

or 

(2)  are  ad  hoc  in  nature  and  fitted  to  reproduce  zero 
temperature  and  classical  results. 
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This  dissertation  has  also  employed  the  dielectric  formalism,  but 
the  electron-electron  interactions  are  incorporated  by  a 
non-perturbati ve  local  field  correction  factor  which  is  obtained  by  a 
quantal  extension  of  classical  fluid  integral  equations.  The 
formalism  asymptotically  reduces  to  standard  statistical  mechanical 
approaches  in  the  classical  limit,  while  in  the  completely  degenerate 
limit  it  reduces  to  a  nonlocal  exchange-correlation  functional 
approximation  within  the  Kohn-Sham  (density  functional)  theory.  The 
validity  of  this  approximation  is  reflected  in  the  good  agreement  we 
have  obtained  at  zero  temperature  with  results  obtained  from  the 

■3Q_42 

Green's  Function  Monte  Carlo  method. 

The  various  approaches  can  be  contrasted  on  several  points. 
First  there  is  a  fundamental  difference  between  those  theories 
utilizing  the  dielectric  formalism  (including  that  of  the  present 
dissertation)  and  the  nonzero-temperature  variational  principle 
utilized  by  Pokrant.    The  dielectric  formalism  must  use  a 
coupling-constant  integration  to  obtain  values  for  the  excess  kinetic 
energy.    This  is  usually  done  implicitly  when  computing  the  free 
energy.    In  the  method  of  Pokrant  the  excess  kinetic  energy  can  be 
computed  explicitly.    Although  it  is  not  given  accurately,  it  is  a 
relatively  unimportant  contribution  to  the  internal  energy,  and  does 
not  affect  values  of  the  free  energy  computed  by  integrating  the 
internal  energy  over  inverse  temperature  at  constant  density.  When 
computing  free  energies  directly,  the  cruder  approximations  in  the 
dielectric  formalism  have  a  machine  time  advantage  over  the  method 
presented  in  this  dissertation.    This  is  because  an  additional  set  of 


-221- 


equations  (the  Kohn-Sham  equations)  appears  within  the 

self-consistency  calculation. 

This  relative  speed  of  numerical  calculation  is  responsible  for 

the  fact  that  the  free  energy  obtained  from  the  cruder  dielectric 

methods  have  been  calculated  over  wide  regions  of  temperature  and 

density.    These  have  been  presented  in  Pade  approximate  form,  where 

the  sole  dependance  upon  the  classical  coupling  parameter,  "gamma,"  is 

modified  by  temperature  dependent  coefficients.    This  method  has  been 

criticized  for  giving  inaccurate  thermodynamic  derivative  quantities; 

this  may  be  attributed  to  parameterizing  the  coupling  in  terms  of 

gamma,  which  is  clearly  inappropriate.    (At  zero  temperature  gamma 

diverges,  while  it  is  known  that  the  free  energy  depends  solely  on  the 
43 

parameter  r^.)      A  more  apt  description  would  be  in  terms  of  the 
parameter 


*3  N 


where 

p  =  (^  ira^)"^ 

^        P       (2vy  2>m 

which  reduces  to  gamma  in  the  classical  limit,  and  is  proportional  to 
r^  in  the  fully  degenerate  limit.    Then  the  free  energy  would  be 
independent  of  a  second  parameter  (needed  to  specify  the  temperature 
and  density)  in  both  extremes,  and  by  inference,  only  loosely 
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dependent  on  it  in  the  intermediate  degeneracy  regime.    Also  this 
description  would  be  more  amendable  to  an  analytic  fit. 

A  second  point  of  contrast  concerns  the  peaked  structure  of  the 
interaction  energy  per  particle  near  absolute  zero.    This  feature  is 
not  present  in  the  method  of  Pokrant.    The  relevant  data  are  not 
directly  reported  in  the  calculations  of  Tanaka  et  al .  or  Dandrea  et 
al .    However  Dandrea  et  al .  do  report  a  feature  that  we  found  to  be 
directly  associated  with  the  hump.    Starting  at  absolute  zero  and 
increasing  temperature,  we  find  that  g(r  =  0),  after  initially 
decreasing  to  a  minimum  value  (located  at  the  temperature 
corresponding  to  the  maximum  of  the  interaction  energy),  then 
increases  for  more  weakly  coupled  plasmas. 

Thirdly,  the  dielectric  method  presented  in  this  dissertation 
admits  modifications  of  the  form  of  the  fluid  equation  obtained  in  the 
classical  limit.    In  classical  statistical  mechanics,  a  blending  of 
fluid  equations  was  found  to  be  an  efficient  algorithm  for  satisfying 
the  compressibility  sum  rule.    Satisfying  this  constraint  generally 
improved  the  accuracy  of  the  thermodynamic  quantities  with  respect  to 
experimental  values.    In  this  research  we  found  that  the  utility  of 
such  a  procedure  was  limited  to  the  near  classical  regime,  but  is 
totally  absent  in  other  dielectric  formalisms. 

Lastly  it  should  be  noted  that  the  method  of  Pokrant  allows  a 
reformulation  of  the  electron  gas  into  a  multispecies  description  of 
spin-up  and  spin-down  particles.    It  has  been  shown^^^  that  in 
statistical  mechanical  descriptions  of  real  plasmas  (i.e.  consisting 
of  positive  ions  and  electrons)  the  resolution  of  the  constituents 
into  spin-  differentiated  species,  while  not  appreciably  affecting  the 


-223- 


charge  structure  factors,  modifies  the  calculated  internal  energies, 
chemical  potentials,  and  compressibilities.    The  extension  of  the 
dielectric  formalism  presented  in  this  dissertation  to  resolve  spin- 
differentiated  species  remains  an  open  avenue  of  research. 
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